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ABSTRACT 


This  thesis  deals  with  the  problem  of  finding  the 
minimum  cost  path  between  two  points  on  a  plane,  for  which 
an  inhomogeneous  cost  is  defined.  It  is  a  problem  in  the 
calculus  of  variations  and  requires  the  minimization  of  a 
cost -di stance  integral.  The  solving  of  the  resulting  Euler- 
Lagrange  equations  is  avoided  by  representing  the  integral 
with  an  objective  function  which  can  then  be  minimized  using 
a  variety  of  mathematical  optimization  techniques.  The 
decision  variables  for  this  objective  function  define  the 
location  of  the  links  of  a  descritized  path,  and  therefore  a 
minimum  of  the  function  provides  the  minimum  cost  path.  This 
method  can  be  used  for  different  cost  surfaces  and  for  paths 
subjected  to  both  linear  and  non-linear  constraints. 

The  case  studies  presented  consider  a  number  of 
surfaces  and  constraints.  The  cost  surface  represented  by 
C ( x , y ) =Kexp ( - Ax 2 -By 2 )  is  used  as  a  theoretical  model  to 
demonstrate  the  interaction  between  cost  and  distance. 
Minimum  time  paths  for  the  velocity  surfaces  of  Wardrop 
(1969),  and  Angel  and  Hyman  (1972b)  are  calculated.  In 
addition,  three  cases  deal  with  imposing  linear  and  non¬ 
linear  constraints  on  the  path. 

Current  literature  explores  the  use  of  transformation 
methods  to  derive  minimum  cost  paths.  Puu  (1978b)  considers 
transforming  the  cost  surface  into  another  surface  which 
upon  distance  equals  cost.  The  minimum  distance  path  defined 


. 


as  the  geodesic  is  shown  to  represent  the  minimum  cost  path. 
The  relationship  between  the  geodesic  and  the  minimum  cost 
path  is  explored  and  discussion  of  the  transformation  method 
is  provided. 

It  is  generally  concluded  that  solving  path  related 
problems  by  minimizing  objective  functions  is  realistic  and 
feasible. 
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Chapter  1  The  Minimum  Cost  Path  Problem 


1 . 1  Introduct ion 

Locating  optimal  paths  in  geographic  space  is  a  problem 
often  encountered  by  man  and  much  energy  is  directed  towards 
its  solution.  Nearly  every  person  undertaking  travel  seeks 
to  optimize  some  factor  such  as  travel  time,  distance,  cost, 
or  even  available  scenery.  Attempts  at  path  optimization  are 
evident  in  the  type  and  location  of  transportation 
facilities  man  has  constructed.  From  the  short  cuts  seen  in 
student  travel  patterns  on  university  campuses,  to 
interconnected  modern  airline  routes,  man  attempts  to  find 
more  efficient  tr anspor tat i on  routes.  With  increasing  energy 
costs  of  providing  both  facilities  and  transport,  the 
optimal  location  of  transpor tat  ion  routes  is  of  increasing 
interest  to  the  geographer,  planner,  and  engineer. 

Optimal  paths  have  been  observed  in  the  movement  of 
nearly  all  natural  phenomena  such  as  the  passage  of  light 
through  various  media,  the  meandering  of  streams,  animal 
paths  of  least  energy  expenditure,  and  avalanche  paths  on 
mountain  slopes.  In  the  past  the  theoretical  analysis  of 
these  paths  has  been  left  more  to  the  physicist  and 
mathematician  than  to  the  geographer.  The  mathematical 
analysis  required  to  explain  the  location  of  optimal  paths 
has  been  beyond  the  interest  and  capabilities  of  most 
geographers.  Current  geography,  however,  requires  the 
methods  of  other  disciplines  to  advance  its  own  theoretical 
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base.  The  main  body  of  theory  concerning  minimum  cost  paths 
occurs  in  physics  and  mathematics.  As  shown  by  Angel  and 
Hyman  (1976)  and  more  recently  by  Puu  (1978  a,b,c),  the 
application  of  optimal  path  theory  is  of  increasing  concern 
to  geographers. 

This  study  is  concerned  with  the  methodology  of  finding 
optimal  paths  connecting  two  points  on  a  plane.  Optimal,  in 
this  sense,  is  defined  in  terms  of  minimizing  some  function 
of  transportation  cost.  The  design  of  this  thesis  is  as 
follows.  Chapter  1  describes  the  minimum  path  problem. 
Chapter  2  develops  the  theoretical  framework  for  finding 
minimum  cost  paths  and  provides  a  brief  review  of  the 
literature.  Chapter  3  explains  the  techniques  of 
mathematical  optimization  used  to  find  the  minimum  cost 
path.  Examples  of  both  unconstrained  and  constrained  minimum 
cost  paths  are  presented  in  Chapter  4.  The  relationship 
between  geodesics  and  minimum  cost  paths  is  examined  in 
Chapter  5.  The  concluding  Chapter  6  provides  a  discussion  on 
the  further  development  of  minimum  cost  path  methodologies. 
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1.2  The  minimum  cost  path 

The  optimal  path  problem  concerns  locating  a  path  of 
minimum  cost  connecting  two  points  on  a  plane.  Before  such  a 
path  can  be  calculated  two  criteria  must  be  met.  First,  the 
path  must  be  free  to  locate  anywhere  on  the  Euclidean  plane 
representing  a  general  physical  surface.  Secondly,  a  cost 
surface  which  provides  a  transportation  cost  for  every  point 
on  the  plane  must  be  defined.  Many  types  of  transpor tat  ion 
cost  variables  have  been  modelled  and  these  can  represent 
the  cost  surface.  Angel  and  Hyman  (1976)  use  travel  velocity 
as  the  inverse  of  a  cost  surface  variable,  Turner  (1971) 
uses  a  combination  of  construction  and  land  acquisition  cost 
variables,  and  Warntz  (1965)  uses  population  potential  as  a 
land  acquistion  cost  variable.  Conceivably,  a  combination  of 
travel,  construction  and  time  costs,  such  as  Wilson's  (1974) 
"generalized  cost  function",  could  be  incorporated  into  the 
cost  surface.  This  research  does  not  concentrate  on  the 
development  of  cost  surfaces;  it  assumes  them  to  be  given. 
Rather,  it  pursues  the  problem  of  finding  minimum  cost  paths 
for  such  surfaces. 

Once  a  cost  surface  and  a  path  have  been  defined,  the 
total  cost  of  the  path  can  be  calculated.  For  the 
homogeneous  or  constant  valued  cost  surface  the  total  cost 
is  calculated  by  multiplying  the  distance  of  the  path  by  the 
cost  per  unit  distance.  The  cost  of  a  path  for  the 
inhomogeneous  or  non  constant  valued  cost  surface  is 
calculated  similarly.  Rather  than  multiplying  the  total 
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distance  of  the  path  by  one  cost,  each  part  of  the  path  must 
be  multiplied  by  the  cost  attributed  to  that  part  and  summed 
for  the  total.  This  calculation  is  equivalent  to  finding  the 
perpendicular  cross-sectional  area  above  the  path  which 
extends  between  the  plane  and  the  cost  surface. 

On  a  geographic  plane  for  which  a  homogeneous  cost 
surface  is  defined  the  minimum  cost  paths  are  straight 
lines.  Figure  1  demonstrates  two  paths  for  a  homogeneous 
cost  surface.  The  constant  value  of  this  surface  is  given  by 
the  formulation  C(x,y)=K,  where  x  and  y  are  coordinates  and 
K  is  a  constant  equal  to  5.0.  The  parameters  "A",  "B" ,  and 
"BC"  which  are  all  equal  to  0.0  are  used  to  formulate  more 
complex  cost  surfaces  discussed  later.  Two  paths  are 
demonstrated  in  this  figure.  Path  1  consists  of  two  links 
with  a  total  length  "DIST"  of  7.80  and  a  total  cost  "COST" 
of  39.02.  Path  2  connects  the  two  points  with  a  straight 
line  of  length  5.66  and  cost  of  28.28.  Each  of  these  costs 
can  be  derived  by  multiplying  the  distance  by  the  constant 
value  of  the  cost  surface.  Figure  2  illustrates  the  area  or 
total  cost  of  each  path.  Path  2,  the  straight  line,  is  the 
path  with  the  least  area  and  cost. 

Minimum  cost  paths  for  inhomogeneous  cost  surfaces  are 
usually  not  straight.  Figure  3  shows  two  paths  on  an 
inhomogeneous  cost  surface.  The  surface  is  projected  onto 
the  plane  by  an  isomap  in  which  isocost  lines  join  points  of 
equal  cost  per  unit  distance.  In  addition  to  the  variable 
cost  surface  each  path  is  subjected  to  a  constant  base  cost 
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Figure  1  Two  paths  for  a  homogeneous  cost  surface 
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point  for  the  paths  in  Figure  1 
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Figure  4  Profile  of  cost  versus  distance  from  one  end  point 
for  paths  in  Figure  3 
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shown  by  the  symbol  "BC"  in  Figure  3.  This  base  cost 
represents  cost  which  can  be  applied  as  a  constant 
throughout  the  length  of  the  path,  such  as  those  arising 
from  road  surfacing,  the  cost  of  wire,  or  even  energy  usage. 
The  isocost  lines  do  not  reflect  the  base  cost.  The  total 
cost  (area)  of  a  path  is  derived  from  two  costs;  the  base 
cost  multiplied  by  the  path  distance  and  the  cost 
contributed  by  the  cost  surface.  For  instance,  if  costs  are 
in  dollars  and  distance  is  in  Kilometers,  the  direct  path  1 
is  5.66  Km.  long  and  costs  $30.60  at  an  average  cost  of 
$5.41  per  Km.  Path  2,  which  avoids  the  high  cost  center,  is 
10.42  Km.  long  with  a  total  cost  of  $12.25  giving  an  average 
cost  of  $1.18  per  Km.  In  the  cost  versus  distance  plot  in 
Figure  4  the  area  for  path  2  is  less  than  that  for  path  1 
demonstrat i ng  that  path  2  is  less  costly  than  path  1.  Even 
though  path  2  is  of  greater  distance  than  path  1,  the  lower 
cost  per  unit  distance  results  in  a  lower  total  cost. 

Intuitively,  minimum  cost  paths  follow  lines  of  least 
resistance.  The  inclination  for  people  to  taKe  short  cuts 
supports  the  least  resistance  principle.  Even  if  the  cost 
criterion  is  changed  from  distance  or  travel  time  to  one  of 
comfort,  the  least  resistance  principle  still  holds. 
Mountainous  landscapes  result  in  roads  that  wind  through 
mountain  passes.  The  slope  and  sinuosity  of  such  roads  is  a 
measure  of  the  resistance  provided  by  the  landscape. 
According  to  Burghardt ( 1 969 )  initial  primitive  roads  located 
to  i nter -connect  administrative  centers  in  the  most 
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efficient  manner  possible.  The  Banff-Jasper  Highway 
connecting  the  two  towns  was  a  dirt  road  constrained  to 
locations  that  allowed  for  easy  river  crossings  and 
avoidance  of  steep  gradients.  The  location,  however,  was 
mostly  through  trees,  avoiding  the  more  scenic  locations. 

The  upgraded  highway,  with  a  higher  construction  budget, 
provided  a  more  scenic  route.  If  scenery  were  considered  to 
be  the  optimization  criterion  the  new  highway  could  be 
considered  to  be  more  cost  effective.  Optimum  paths  within 
urban  environments  often  relate  costs  to  travel  time.  Ring 
roads  around  urban  centers  are  located  to  trade  off  the  cost 
of  land  acquisition  with  the  ability  of  the  road  to  ease 
local  traffic  congestion.  Regardless  of  the  type,  many  paths 
result  from  attempts  to  optimize  some  real  or  perceived 
cost . 


Locating  a  minimum  cost  path  is  not  a  trivial  problem. 
Intuitive  methods  can  result  in  non-optimal  paths,  because 
the  chosen  path  often  accumulates  too  much  cost  by  being 
either  too  long  or  too  short.  In  addition,  optimal  paths 
often  trade  distance  and  cost.  Sometimes  paths  traversing  a 
short  high  cost  area  are  more  economical  than  those  taking  a 
longer  but  less  costly  route.  The  reverse  may  also  be  true. 
Only  since  Werner  and  Boukidis  (1963)  have  intuitive  methods 
been  given  theoretical  and  mathematical  expression.  The 
provision  of  such  expressions  has  not  reduced  the  problem  to 


a  trivial  one,  for  their  analysis  requires  the  solution  of 
difficult  partial  differential  equations.  Modern  numerical 


methods  and  computers  have  allowed  workers  such  as  Turner 
(1971)  to  analyse  applied  engineering  problems  and  provide 
the  geographer  with  the  tools  to  analyse  a  wide  range  of 
optimal  path  problems. 


Chapter  2  Known  Solution  Methods 


2.1  Early  mathematical  development  of  path  problems 

In  their  fundamental  dissertation  What  i s  Mathemat i cs? , 
Courant  and  Robbins  (1941)  demonstrated  that  early 
mathematics  was  developed  from  the  geometry  of  points, 
lines,  and  surfaces.  In  geometry  many  maximum  and  minimum 
problems  concerning  lines  (paths)  have  been  formulated  and 
the  solution  of  these  problems  has  led  to  the  mathematics 
used  today.  The  following  discussion  presents  some  of  the 
classical  path  problems  Known  for  the  last  twenty  centuries. 

The  Greeks  accepted  without  proof  that  the  shortest 
path  between  two  points  on  a  plane  is  a  straight  line. 
Although  many  theorems  in  the  Euclidean  geometry  depended  on 
this  observation,  the  Greeks  found  no  reason  to  doubt  this 
axiom.  Perhaps  the  Pythagorean  theorem  may  have  served  as  an 
intuitive  proof.  The  Euclidean  geometry  was  not  only  an 
academic  exercise  for  the  Greeks,  but  also  a  tool  for 
explaining  natural  phenomena.  Heron,  an  Alexandrian  living 
in  the  first  century,  showed  that  the  angle  of  incidence  of 
a  reflected  light  ray  was  equal  to  the  reflective  angle. 
Courant  and  Robbins  considered  this  discovery  to  be  "the 
germ  of  the  theory  of  geometrical  optics". 

The  contribution  towards  minimum  path  theories  by 
problems  concerning  light  rays  is  extensive.  During  the  late 
rennai ssance ,  Snell  (1591-1626)  formulated  the  law  of 
refraction  which  describes  the  path  location  of  a  light  ray 
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as  it  passes  from  one  medium  to  another.  His  law  states  that 
sin  A/sin  B=u/v,  where  A  and  B  are  the  angles  the  light  ray 
makes  with  the  perpendicular  interface,  and  u  and  v  are  the 
velocities  of  the  light  in  the  two  mediums.  Later  in  the 
seventeenth  century  Fermat  used  the  calculus  to  prove  that 
the  refracted  path  of  Snell  was  indeed  the  shortest  time 
path  through  the  two  media.  Fermat  extended  this  analysis  to 
light  traveling  through  many  media  of  infinitely  small 
thicknesses  and  developed  his  "principle  of  geometric 
optics".  Losch  (1954)  used  this  refraction  law  to  describe 
the  location  of  the  least  cost  transportation  route  in  a 
sea/land  interface. 

Fermat  was  not  the  only  mathematician  who  contributed 
to  the  calculus  and  minimum  path  problems  in  the  17' th 
century.  The  famous  "brachi stochrone"  problem  posed  by 
Johann  Bernoulli  for  his  " i ncompetant "  older  brother  John 
created  even  more  interest  in  minimum  path  problems. 

"Imagine  a  particle  constrained  to  slide  without  friction 
along  a  certain  curve  joining  a  point  A  to  a  lower  point  B. 
If  the  particle  is  allowed  to  fall  under  influence  of 
gravity  alone,  along  which  such  curve  will  the  time  required 
for  the  descent  be  least?"  (Courant  and  Robbins,  page  379 
quoting  ACTA  ERUDITRUM  a  17'  th  century  scientific  journal). 

The  problems  of  Fermat  and  Bernoulli  were  recognized  in 
their  time  as  not  being  solvable  by  differential  calculus, 
which  required  known  functions  (or  curves)  for  its 
operation.  Euler  and  Lagrange  (1736-1813)  found  a  general 
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method  for  determining  paths  which  minimize  some  variable(s) 
such  as  time  in  the  "brachi stochrone"  problem  or  distance  in 
Fermat's  principle.  Their  method,  called  the  calculus  of 
variations,  was  one  of  the  most  important  developments  for 
the  applied  mathematics  and  can  be  found  in  the  solution 
methods  for  many  optical,  mechanical,  physical,  and 
cartographic  problems. 

Geodesics  have  been  in  use  in  cartography  for  many 
centuries,  yet  their  contribution  to  spatial  analysis  has 
until  recently  been  limited  mainly  because  the  calculation 
of  geodesics  on  surfaces  other  than  a  plane  or  a  sphere 
requires  a  calculus  of  variations  approach.  The  geodesic  is 
the  shortest  path  between  two  points  on  a  surface.  Geodesics 
form  the  basis  for  understanding  and  defining  many  non- 
Euclidean  geometries.  Bernard  Riemann  (1826-1866)  discovered 
the  elliptic  geometry,  in  which  the  geodesics  are  straight 
lines.  The  earth  is  an  ellipsoid  and  the  great  circle  curves 
or  geodesics  on  the  surface  are  defined  as  straight  lines  in 
the  elliptical  geometry. 


2.2  General  mathematical  formulation 

Existing  general  mathematical  formulations  for  defining 
minimum  cost  paths  may  be  separated  into  two  types:  first, 
continuous  path  formulations  requiring  continuous  cost 
surfaces  and  secondly,  discrete  path  formulations  which  have 
been  developed  for  both  discrete  and  continuous  cost 
surfaces.  Werner  and  Boukidis  (1963)  provide  the  most 
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general  formulation.  Their  development  is  incorporated  into 
the  following  discussion. 


2.2.1  Continuous  path  formulations 

For  a  2-dimensional  geographic  plane  described  by  x  and 
y  coordinates,  the  cost  at  any  point  p(x,y)  is  given  as: 

C  =  f(x,y)  .  .  .  (  1  ).  Def i ne  a  continuous  path  s  between  two 
points  on  the  geographic  plane  as  an  infinite  number  of 
connected  links  of  length  ds ,  where 

ds  =  ^(dx)2  +  (dy)2  . (2) 

The  total  cost  T  of  a  path  constructed  between  two  points 
Po  (x0  ’V0  )  and  pf  *xf  »yf  )  is  given  by: 


T(s) 


r p< 

=/C(x, 


y  )ds , 


=/c('x,  y)^ 


1+ (dy/dx ) 2  dx 


(3) 


Werner  and  Boukidis  (1963)  recognized  that  the  minimization 
of  (3)  requires  a  calculus  of  variations  approach.  Equation 
(3)  can  be  generalized  as: 


Minimize 


.y.y 


/ 


)dx , 


(4) 


where  y'  =dy/dx .  The  solution  to  this  equation  is  given  by 
the  Euler  partial  differential  equation 

c»F/£y  -  d(^F/ay'  ) /dx  =  0 . (5) 

Werner  and  Boukidis  (1963)  provide  a  detailed  development. 
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A  more  rigorous  solution  description  for  the  general 
equation  (4)  was  given  by  Howard  et  aj_.  (1968).  The 
coordinate  variables  are  placed  in  parametric  form  as  x(t) 
and  y(t).  The  parameter  t  could  be  defined  as  travel  time. 
Then,  equation  (4)  becomes: 


/< 

Minimize  T=  / F ( x ( t ) , y ( t ) , x ( t ) , y ( t ) ) dt , 


(6) 


where  x(t)=dx/dt,  and  y(t)=dy/dt.  The  solution  of  (6)  is 
given  by  the  Euler-Lagrange  equations: 


xC(  x  ,y )  \  =^x2+y2  d_F 

Jx2+y2  /  5x . (7.1) 

yC ( x , y )\  =^*2+y2  3F 

fi2^2  )  sy . (7.2) 

which  are  the  conditions  for  an  optimal  trajectory  as 
developed  by  Pontryagin  (1961).  This  approach,  called  the 
optimum  curvature  principle  by  Howard  et  aj.  (1968),  was 
developed  for  the  optimum  location  of  highways. 


2.2.2  Discrete  formulations 

Generally  not  all  cost  surfaces  or  paths  can  be 
represented  in  continous  form.  Werner  and  Boukidis  (1963) 
suggest  that  continuous  problems  are  specific  cases  of  the 
more  general  discrete  formulation.  The  integral  of  equation 
(3)  becomes  a  summation.  The  total  cost  from  point  P(x0,y0  ) 
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to  P ( x  . y  )=P(x  y )  is  defined  as: 
n+|  n+1  f  f 

n  - 1 

T=^C  (  x  ,  y )  J  Ax.  2  +  Ay.2  . (8) 

i  =  0 

where  n  =number  of  links,  Ax.  =  x. ,  ,-x.  ,  and  Ay.  =y  -y  .  As 

i  '+1  '  J\  i+i  J\ 

Ax.  ,  Ay.  ,  approaches  zero  (8)  approaches  the  continuous 
formulation  (3).  The  path  s  is  defined  by  the  points 
P ( x  , y.  ) ,  i  =  1,...,n.  This  discrete  for mu  1  at  ion  is  more 
general  than  the  continuous  one  because  of  the  different 
ways  the  cost  surface  can  be  defined.  The  cost  C(x,y) 
defined  at  a  point  P(x,y)  can  be  represented  in  a  number  of 
ways : 

-  by  a  continuous  function  as  defined  by  (1) 

-  by  contiguous  areas  of  constant  cost  (cost  planes) 

-  by  a  matrix  of  costs  representing  values  at  points 
P(x,y)  in  the  region  of  interest.  The  solution  methodology 
appropriate  to  minimizing  (8)  depends  on  the  characteristics 
of  the  cost  surface. 

If  the  cost  surface  C(x,y)  is  smooth  the  Euler-Lagrange 
equations  (7.1)  and  (7.2)  may  be  used  to  find  the  minimum 
path.  Not  all  cost  surfaces  provide  for  simple  analytical 
solutions  to  the  Euler-Lagrange  equations,  so  numerical 
techniques  as  suggested  by  Howard  et  aj_.(1968)  must  be  used. 
Rankin  (1979)  developed  an  "artillary"  method  to  solve  the 
Euler-Lagrange  equations.  His  method  utilizes  the  boundary 
conditions  that  require  the  solution  to  pass  through  both 
end  points.  By  shooting  from  one  point  to  the  other  with 
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different  angles  and  path  lengths,  and  integrating  the 
Euler-Lagrange  equations  using  a  4' th  order  Runga-Kutta 
algorithm,  Rankin's  method  converges  to  an  extremum  path. 
This  path  is  then  examined  for  a  possible  minimum. 

A  general  method  for  finding  minimum  cost  paths  on  cost 
surfaces  defined  by  contiguous  areas  of  constant  cost  has 
been  presented  by  Werner  (1968).  He  extends  Snell's 
refraction  principle  and  the  application  of  Losch  to  a  path 
traversing  a  number  of  cost  planes  and  proves  that  a  unique 
set  of  refraction  angles  exist  for  the  minimum  cost  path. 
Wardrop  (1969)  also  uses  the  refraction  method  in 
calculating  optimum  paths  for  continuous  cost  media. 

Minimum  path  algorithms  provide  solution  paths  for  cost 
surfaces  represented  at  discrete  points.  Each  link  between 
two  points  on  the  geographic  plane  is  given  a  value 
proportional  to  the  cost  of  connecting  the  two  points.  A 
matrix  of  costs  is  constructed,  where  a  node  (or  element) 
(i,j)  of  the  matrix  provides  the  cost  of  the  link  connecting 
the  point  P(x.  ,y.  )  with  P  ( ,y^  ).  Two  types  of  algorithm  are 
described  by  Steenbrink;  the  tree  building  algorithms,  which 
build  the  shortest  tree  from  each  node  to  all  other  nodes, 
and  matrix  algorithms,  which  calculate  the  shortest  paths 
between  all  nodes  and  all  other  nodes  and  store  the  paths  in 
matrix  form.  Moore  (1959)  developed  the  earlier  tree 
building  algorithms.  Following  this,  "once  through"  tree 
building  algorithms  were  published  by  Dijkstra  (1959)  and 
Whiting  and  Hillier  (1960).  Efficient  matrix  algorithms  were 
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developed  by  Floyd  (1962),  Dantzig  (1966)  and  Farbey  et  a  1 . 
(1967).  One  algorithm  developed  by  Murchland  (1967)  allows 
for  a  change  in  one  link  without  recomputing  all  the 
shortest  paths.  Steenbrink  discussed  heuristic  procedures 
which  provide  very  good  solutions  with  efficient  use  of 
computer  resources,  computer  resources.  One  such  algorithm 
developed  by  Murchland  (1967)  allows  for  a  change  in  one 
link  without  recomputing  all  the  shortest  paths.  Steenbrink 
also  mentions  the  use  of  dynamic  programming  as  is 
incorporated  into  the  solution  procedures  described  by 
Turner  (1971).  Goodchi Id ( 1 977 )  recently  developed  heuristic 
procedures  to  find  minimum  cost  paths  through  lattices 
defined  in  an  urban  setting. 

2.2.3  Alternative  Methods 

Alternative  methods  for  finding  optimum  paths  have  been 
developed  by  geographers  and  engineers  in  order  to  solve 
specific  problems  within  their  respective  disciplines. 
Surface  transformation  methods  are  the  most  important  of  the 
alternate  methods  for  geographers,  because  they  have  a  solid 
basis  within  cartography.  The  engineer,  who  is  more 
interested  in  practical  applications  of  the  optimum  path 
methodology,  has  developed  complete  computer  systems  for  the 
analysis  of  cost  data  and  feasible  paths. 

The  soap  film  experiments  of  the  Belgium  physicist 
Plateau  (1801-1883)  provided  general  insights  in  solving 
calculus  of  variations  problems.  Plateau  observed  that  soap 
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films  suspended  between  surfaces  and  or  edges  minimized 
their  areas.  Courant  and  Robbins  (1941)  demonstrate  a  wide 
range  of  geometrical,  topological  and  variational  problems 
solved  by  soap  film  analogies.  Steiner's  problem  and  it's 
three-dimensional  extensions  are  easily  solved.  An 
experiment  to  solve  the  problem  of  this  thesis  could  also  be 
designed  using  soap  films.  A  film  is  suspended  between  the 
closed  area  bounded  by  the  plane,  the  cost  surface,  and  two 
perpendicular  edges  above  the  end  points.  After  the  soap 
film  is  stabilized,  the  contact  by  the  film  on  the  plane 
should  trace  the  minimum  cost  path. 

Minimum  paths  for  radially  symmetric  cost  surfaces  have 
received  recent  attention  by  Angel  and  Hyman  (1972a),  and 
Puu  (1978  a,b,c).  They  have  transformed  the  velocity  (cost) 
surface  into  another  surface  where  upon  the  geodesic 
corresponds  to  the  minimum  time  (cost)  path.  Their  methods 
are  discussed  in  Chapter  5. 

A  systems  approach  employing  a  wide  range  of  algorithms 
has  been  used  for  optimal  route  location  by  engineers  in 
both  North  America  and  Europe.  Howard  e_t  aj_.  (1968)  suggest 
that  their  optimum  curvature  principle  can  be  incorporated 
within  an  engineering  system  such  as  "TIES"  (Total 
Integrated  Engineering  System)  described  by  Schureman  (1965) 
or  the  Massachusetts  Institue  of  Technology  ICES  (Integrated 
Civil  Engineering  System)  as  described  by  Ross  and 
Schemakcer  (1965).  Dynamic  programming  techniques  are  used 
to  solve  the  Euler-Lagrange  equations  in  these  systems.  More 
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recent  systems  such  as  OECD  (Road  Research  Group  in  Paris), 
developed  in  Paris,  and  GCARS  (Generalized  Computer-Aided 
Route  Selection  System)  developed  by  Turner  and  Miles  (1971) 
utilize  a  combination  of  minimum  cost  algorithms  and  dynamic 
programming  methods.  The  GCARS  system  uses  the  once  through 
algorithm  of  Whiting  and  Hillier  (1960)  in  conjunction  with 
an  objective  function  which  consists  of  a  "utility  surface" 
calculated  with  a  combination  of  trend  surfaces  representing 
local  cost  factors  and  the  distance  calculated  from  a 
network  of  interconnected  links. 

2.3  Discussion 

Until  now,  the  existing  methods  have  been  mentioned, 
but  no  attempt  has  been  made  to  evaluate  them  as  to  their 
suitability  to  solve  theoretical  and  empirical  geographic 
problems.  The  following  discussion  considers  the  methodology 
choosen  for  this  thesis.  A  common  set  of  criteria  which  can 
be  used  to  evaluate  these  methods  is  difficult  to  obtain, 
because  the  existing  methods  are  tailored  to  solve  specific 
types  of  problem. 

The  first  of  the  methodologies  is  contained  within  the 
mathematical  development  given  by  Werner  and  Boukidis  (1963) 
and  Howard  et  ad.  (1968).  Although  the  Euler-Lagrange 
equations  (7.1)  and  (7.2)  provide  a  good  theoretical 
definition  of  the  location  of  the  minimum  cost  path,  the 
exact  solution  of  these  equations  is  generally  very 
difficult.  This  is  supported  by  Courant  and  Robbins  (1941). 
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"It  is  usually  very  difficult,  and  sometimes 
impossible,  to  solve  variational  problems  explicitly  in 
terms  of  formulas  or  geometrical  constructions  involving 
simple  elements.  Instead,  one  is  often  satisfied  with 
merely  proving  the  existence  of  a  solution  under  certain 
conditions  and  afterwards  investigating  properties  of  the 
solution.  In  many  cases,  when  such  an  existence  proof 
turns  out  to  be  more  or  less  difficult,  it  is  stimulating 
to  realize  the  mathematical  conditions  of  the  problem  by 
corresponding  physical  devices,  or  rather,  to  consider 
the  mathematical  problem  as  an  i nterpretat ion  of  a 
physical  phenomenon.  The  existence  of  the  physical 
phenomenon  then  represents  the  solution  of  the 
mathematical  problem."  (Courant  and  Robbins,  page  386) 
Howard  et.  a]..  (1968)  provide  a  family  of  solutions  to  the 
Euler-Lagrange  equations  for  an  exponential  cost  model  by 
finding  i ni t i a  1 -va 1 ue  conditions  that  satisfy  the  equations. 
They  suggest  that  solution  conditions  for  the  two-point 
boundary  value  problem  are  very  much  harder  to  find  and  that 
numerical  integration  of  the  objective  function  (6)  is  an 
easier  problem.  Generally,  however,  the  solution  of  these 
equations  for  a  geographic  application  is  too  tedious.  The 
use  of  soap  film  experiments  to  provide  solutions  to  these 
equations  also  has  practical  problems,  because  a  physical 
representation  of  the  cost  surface  is  difficult  to 
construct.  If  constraints  are  placed  on  path  location,  new 
Euler-Lagrange  equations  may  have  to  be  developed  from  first 
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principles.  Rankin  (1979)  mentions  that  considerable 
modi  feat  ion  to  his  method  must  be  made  before  he  can 
incorporate  constraints  into  his  numerical  method.  The  lack 
of  literature  exploring  this  problem  area  indicates  the 
difficulty  with  this  purely  theoretical  approach. 

The  representation  of  both  continuous  and  discrete  cost 
surfaces  by  many  planes  of  constant  cost  values,  as 
suggested  by  Werner  (1968)  demands  the  acceptance  of  two 
assumptions.  Firstly,  Werner  (1968)  assumes  that  the 
geographical  cost  surface  may  be  adequately  represented  by  a 
set  of  these  contiguous  cost  planes.  Secondly,  he  assumes 
that  the  refracted  path  at  the  interface  between  two  cost 
planes  is  acceptable  and  realistic  for  many  types  of  paths. 
For  paths  representing  roads  through  landscapes  with  small 
variations  in  cost  the  method  is  proven  to  be  adequate  by 
Werner  and  Boukidi s ( 1 963 ) .  However,  some  cost  surfaces 
representing  travel  time  within  an  urban  environment,  may  be 
very  difficult  to  approximate  with  constant  cost  planes.  The 
cost  surface  may  have  very  steep  gradients  and  if  the  number 
of  cost  planes  representing  this  surface  is  small,  the 
refraction  angle  of  the  path  may  be  very  large.  By 
increasing  the  number  of  cost  planes  the  path  may  locate  in 
a  different  region.  Werner's  method  also  requires  that  the 
set  of  planes  to  be  traversed  by  the  minimum  cost  path  be 
known  before  the  refraction  angles  can  be  calculated.  This 
requirement  may  be  met  by  the  solution  methods  proposed  in 
this  thesis.  In  addition  Werner  suggests  that  constraints 
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must  be  considered,  but,  he  does  not  incorporate  them  into 
his  solution  methodology.  Werner  concludes  that  his 
particular  "model  provides  only  basic  concepts  on  which  to 
develop  a  theoretical  approach". 

The  most  general  of  all  solution  methodologies  are 
found  in  minimum  path  algorithms.  Once  the  surface  is 
represented  by  point  to  point  costs  a  variety  of  algorithms 
can  be  used  to  calculate  the  minimum  cost  paths.  The  use  of 
these  algorithms  by  engineers  indicates  their  advantage  in 
being  suitable  for  solving  practical  problems.  The 
disadvantages  of  using  minimum  path  algorithms  are 
threefold.  Firstly,  for  very  large  problems  the  computation 
costs  of  some  algorithms  are  unreasonable,  even  though  a 
variety  of  heuristic  techniques  can  be  used  to  provide 
efficiency  (see  Steenbrink  (1974),  Chapter  7).  Secondly,  the 
solution  vector  is  restricted  to  the  grid  points  defining 
the  cost  surface,  hence,  the  minimum  cost  path  is  no  more 
precise  than  the  density  of  the  grid.  The  third  disadvantage 
concerns  constraining  the  path  to  certain  criteria  such  as 
slope,  radius  of  curvature,  or  other  factors.  Either  the 
minimum  path  algorithm  or  the  cost  matrix  must  be  modified 
to  incorporate  this  type  of  constraint. 

The  systems  TIES  and  GCARS  have  been  used  to  resolve 
many  empirical  engineering  problems  concerning  minimum  cost 
paths.  The  GCARS  method  uses  the  minimum  path  algorithms 
discussed  above  and  mentions  the  same  concern  for 
computational  efficiency  with  large  grids.  The  TIES  system 
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uses  minimum  path  algorithms  in  conjunction  with  numerical 
integration  of  the  cost  integral  to  define  alternative 
minimum  cost  paths.  The  expansion  of  these  systems  into  a 
general  framework  which  is  suitable  for  theoretical 
geographic  application  is  a  software  evaluation  problem 
beyond  the  scope  of  this  thesis. 

2.4  The  mathematical  optimization  formulation 

Aside  from  the  GCARS  system  described  by  Turner  (1971), 
there  have  been  very  few  attempts  made  to  use  a  mathematical 
optimization  approach  to  solve  the  basic  minimum  cost  path 
problem.  Steenbrink  (1974)  suggested  that  a  discrete 
formulation  such  as  presented  by  an  equation  similar  to  (3) 
can  be  very  easily  solved  using  dynamic  programming  if  the 
path  can  be  represented  as  a  sequence  of  connected  points. 
The  use  of  mathematical  optimization  to  solve  minimum  cost 
path  problems  especially  related  to  road  location  was 
discussed  in  the  two  symposia  of  the  Planning  and  Transport 
Reasearch  and  Computation  (1969  and  1971)  and  initially  by 
Werner  and  Boukidis  (1963). 

A  mathematical  optimization  formulation  of  the  minimum 
path  is  given  by: 

n- 1 

Minimize  F  )  D . . (9) 

i  =  0 

where  C(Pj  )=C(xs  ,y.  )  and  D-  =^(  xj+1-x.  )  2  +  (y.+1  -y.  ) 2  ,  and  n 
equals  the  number  of  links  connecting  P(xo,yo  )  and  P  ( xf  ,yf  ). 
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Equation  (9)  is  Known  as  an  objective  function,  the  meaning 
of  which  can  be  expressed  as:  "which  points  P  , i=1,...,n 
minimize  the  value  of  F  (or  the  cost)".  A  general  treatment 
of  this  objective  function  is  given  in  the  next  chapter. 

The  use  of  mathematical  optimization  has  merit.  By 
minimizing  the  discrete  version  of  equation  (3),  there  is  no 
need  to  solve  the  Euler-Lagrange  equations  (7.1)  and  (7.2). 
The  objective  function  is  flexible.  The  general  cost  surface 
C(x,y)  may  contain  any  combination  of  cost  surfaces  as  long 
as  each  surface  is  continuously  differentiable  for  every 
point  (x,y).  The  path  defined  by  n  links  can  also  be 
modified  to  include  costs  dependent  on  only  the  path.  Both 
linear  and  nonlinear  constraints  imposed  on  the  path  can  be 
dealt  with  according  to  well  developed  optimization 
methodology.  Finally,  different  objective  functions  can  be 
developed  to  accommodate  a  wide  variety  of  path  minimization 
problems . 


Chapter  3  The  Objective  Function:  Imp  1 ementaton  and  Solution 


3.1  Mathematical  optimization 

The  minimization  of  the  objective  function  (9)  may  be 
achieved  through  the  use  of  non-linear  mathematical 
optimization  techniques.  A  general  introduction  to 
mathematical  optimization  can  be  found  in  AoK i ( 1 9 7 1 ) , 

Steenbr ink( 1974 ) ,  or  in  a  variety  of  text  books.  Steenbrink 
formulates  the  general  mathematical  optimization  problem  as 
fol lows : 

"It  is  desired  to  determine  values  for  n  variables 

x1  . xn  in  such  a  way  that  the  value  of  a  function 

of  these  variables  ( F ( x1  ,  .  .  . , xn ) )  is  as  1 arge  or  as 
small  as  possible.  The  variables  x1  ,  .  .  .  ,  xn  are  called 
decision  variables  or  instrument  variables.  The 

function  F ( x1  . xn )  is  called  the  objective 

function.  The  objective  function  must  be  maximized  or 
minimized . 

Moreover,  there  (may)  exist  certain 
relationships  between  the  decision  variables  and/or 
the  decision  variables  or  functions  of  the  decision 
variables  which  must  satisfy  some  inequalities  or 
equalities.  These  are  called  constraints:..." 

The  minimization  problem  is  generally  formulated  as: 

Minimize  F  ( x1  ,  .  .  .  ,  xn  ) . (10) 

subject  to  gj  ( x1 . xn ) <0 ,  j= 1 , nc , 


27 


28 


and  h  (  x1  ,  .  .  .  ,  xn  )  =0  ,  k=  1  ,  neq  ;  . (11) 

where  nc  equals  the  number  of  inequality  constraints  and  neq 
equals  the  number  of  equality  constraints. 

The  problem  is  expressed  in  vector  notation  as: 

minimize  F ( X  )  ,  s.t.  G(X)<0,  H ( X ) =  0 ,  . (12) 

where  X=x1  . xn .  The  conditions  for  minimizing  this 

objective  function  are  seperated  into  unconstrained  and 
constrained  types. 

The  conditions  for  deriving  a  minimum  of  an 
unconstrained  objective  function  are  now  examined.  Define  X* 
as  an  optimal  vector  of  decision  variables  that  minimizes 
(10).  The  necessary  conditions  for  the  minimum  X*  require 
that  every  partial  derivative  of  F  with  respect  to  x  , i = 1 , n 
be  equal  to  zero.  The  vector  of  partial  derivatives  is 
called  the  gradient  and  is  given  by: 

VFUMdF/dx, . ^F/8xn  ) . (13) 

For  a  convex  F,  X*  is  a  local  and  a  global  minimum  if  and 
only  if  the  gradient  is  equal  to  zero.  F  is  convex  if  for 
any  X,  ,X2 

F(XX1+(  1-X)X2  )  <XF(X1  )  +  (  1  -  X )  F  (  X2  )  ,  if  0<X<1. 

If  F  is  strictly  convex  and  the  gradient  is  equal  to  zero, 
then  X *  is  a  unique  global  minimum.  F  is  strictly  convex  if 
for  any  X1 , X2 

F(AX1 +( 1-A)X2 )<XF(X1  )  +  ( 1-X)F(X2 ) ,  if  0<X<1. 

The  sufficient  conditions  for  a  local  minimum  can  be 
tested  by  deriving  the  Hessian  matrix  of  dimension  (nXn). 
This  matrix  is  calculated  by  taking  all  the  partial 
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derivatives  of  the  gradient  at  X*  which  are  given  by: 

HESS  (  X*  )  =  ^  F  2/  ax.  d  Xj  . (14) 

where  i = 1 , n  and  j=1,n.  A  local  minimum  can  be  declared  if 
the  Hessian  matrix  is  positive  definite.  A  symmetric  matrix 
is  positive  definite  if  and  only  if  all  eigenvalues  are 
greater  than  zero.  A  positive  semidefinite  matrix  requires 
all  eigenvalues  to  be  greater  or  equal  to  zero.  If  the 
Hessian  matrix  is  positive  semidefinite  at  X*,  no  decision 
on  whether  X*  is  a  local  minimum  can  be  made.  Instead  the 
eigenvalues  of  the  Hessians  evaluated  in  the  neighborhood  of 
X*  must  be  examined.  If  the  eigenvalues  of  these  Hessians 
are  all  non-negative,  then  a  local  minimum  can  still  be 
claimed.  A  more  complete  discussion  of  the  conditions 
required  for  local  and  global  minimization  is  found  in 
Steenbr i nk ( 1 974 ) ,  AoK i ( 1 97 1 ) ,  and  Hauer(1974). 

For  general  minimization  problems  with  constraints,  the 
objective  function  (10)  and  constraints  (11)  are  combined 
with  the  use  of  Lagrangian  multipliers  to  give  the 
unconstrained  formulation: 

minimize  K(X,L,U)=F(X)+LTG(X)+UTH(X),  . (15) 

where  L=1  ,  .  .  .  ,  lnc  and  U  =  u1  ,  .  .  .uReq  .  Given  that  X*  is  a 
local  minimum  which  satisfies  the  equality  and  inequality 
constraints  given  by  (11)  and  at  X*  a  suitable  constraint 
qualification  holds,  then  there  exist  1 j  > 0 ,  i  =  1 , nc ,  such 
that  the  following  conditions  hold: 


Kx  (  X*  ,  L  ,  U  )  =FX  (X*)  +  LTGx(X*)+UTHx(X*)=0, 
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and  LTG(X*)=0 . (16) 

where  L  and  U  are  as  defined  in  (15).  These  equations  are 
Known  as  the  Kuhn- Tucker  conditions  defined  by  Kuhn  and 
Tucker  (1951). 

Aok i ( 1971 )  classified  the  methods  of  minimizing  or 
maximizing  an  objective  function  into  three  categories: 

(1)  Methods  using  only  the  functional  values,  called 
direct  methods. 

(2)  Methods  making  use  of  the  first-order  derivatives 
as  we  1 1 . 

(3)  Methods  which  also  require  knowledge  of  second- 
order  derivatives.  Aok i ( 1 9 7 1 ) . 

Generally,  the  solution  methods  used  depend  on  the 
characteristics  of  the  objective  function.  Direct  methods 
are  used  for  objective  functions  for  which  the  first  and 
second  derivatives  are  difficult  to  calculate.  First-order 
derivative  methods,  commonly  known  as  gradient  methods, 
offer  the  most  flexibility  for  solving  problems  which  have 
difficult  solutions  and  imposed  constraints.  Methods 
utilizing  second  derivatives  require  an  objective  function 
and  constraints  with  existing  second-derivatives.  Most 
methods  locate  stationary  points  and  the  convexity 
conditions  must  be  examined  at  these  points  before  local  or 
global  minima  can  be  claimed. 
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Several  optimization  techniques  using  gradient  methods 
have  been  programmed  by  Hauer (1974).  These  methods  utilize 
the  first  order  derivatives  of  the  objective  function  in  an 
iterative  procedure.  An  initial  approximate  solution  X  and  a 
search  direction  equal  to  the  negative  reciprocal  of  the 
gradient  is  given  to  the  program.  The  next  approximate 
solution  X  is  found  with  the  recursion 
X ( k+ 1 ) =X ( K ) +t ( K ) D ( k ) ,  K=0 , 1 , 2 , 3 , . . . . ;  where  K  is  the 
iteration,  D  is  the  feasible  search  direction  and  t  is  a 
scalar  determining  the  step  size.  After  each  iteration  a  new 
feasible  search  direction  is  determined  from  the  gradient  at 
X(k).  The  procedure  stops  at  a  stationary  point  if  the 
following  conditions  occur: 

(1)  if  the  gradient  beomes  sufficiently  close  to  zero 
for  X ( k )  in  an  unconstrained  region 

(2)  if  any  search  direction  results  in  a  non-decrease 
in  the  value  of  the  objective  function 

(3)  if  no  significant  decrease  in  F  can  be  achieved 
in  a  fixed  number  of  iterations. 

For  constrained  minimization  the  Kuhn-Tucker  conditions  at 
the  stationary  point  are  examined.  If  the  conditions  (16) 
are  satisfied,  and  if  the  objective  function  and  the 
inequality  constraints  are  convex,  while  the  equality 
constraints  are  affine,  the  stationary  point  is  declared  a 
global  minimum  at  the  constraints.  For  unconstrained 
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minimization,  the  gradient  and  the  Hessian  at  the  stationary 
point  are  examined  for  a  local  minimum. 

The  advantages  of  using  the  gradient  methods  for  our 
application  can  be  summarized  as  follows: 

(1)  The  objective  function  is  smooth  and  convergence 
occurs  in  a  few  iterations. 

(2)  The  computer  algorithms  available  concentrate  on 
gradient  methods. 

(3)  Constraints  can  be  incorporated  into  the  solution 
procedures . 

A  more  complete  description  of  these  methods  is  found  in 
Aoki(1971)  and  Hauer(1974). 

3.2  Modelling  the  minimum  cost  path  as  an  objective  function 

n-1 

The  objective  function  (9)  F=  ^C(P.  )D.  may  clearly 

i  =  0  1 

represent  the  theoretical  problem,  but  it  is  not  a 
particulary  useful  function  for  optimization.  The  locational 
variables  P1  ,  .  .  .  ,  Pn  as  defined  in  section  2.4  are  too 
general  and  must  be  represented  by  a  set  of  decision 
variables  which  can  be  manipulated.  There  are  many  possible 
ways  of  combining  the  x  and  y  coordinates  of  the  locational 
variables  into  this  set  and  not  all  methods  provide  good 
solutions  of  the  objective  function.  This  section  addresses 
the  problem  of  choosing  a  set  of  decision  variables  which 
provide  optimum  paths  and  at  the  same  time  satisfy  the 
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operational  requirements  of  the  solution  methodology. 

3.2.1  The  two  link  path  and  the  objective  function 

The  technique  of  isodapanes  developed  by  Weber  (1909) 
may  be  used  to  demonstrate  the  objective  function. 
Isodapanes  are  lines  of  equal  total  cost  from  two  fixed 
points  on  the  plane  to  a  variable  third  point.  If  the 
objective  function  (9)  is  restricted  to  a  single  locational 
variable,  then  the  contours  of  the  function  correspond  to 
Weber's  isodapanes.  Isodapanes  were  used  by  Smith  (1966)  to 
find  locations  which  minimize  the  total  transportation  cost 
of  raw  materials  and  finished  products.  Similary,  the 
contours  can  be  used  as  a  graphical  tool  to  find  the 
location  on  the  plane,  where  the  objective  function  is  at  a 
minimum.  The  general  form  of  the  objective  function  with  a 
single  locational  variable  P(x,y)  is  given  by: 


(17) 


where  Co(x,y)  is  the  average  cost  above  the  first  link  and 
Cf(x,y)  is  the  average  cost  above  the  last  link  in  the  path. 
This  average  cost  is  explained  in  the  Section  3.2.2.  The 
decision  variables  of  this  objective  function  are  the  x  and 
y  coordinates  of  the  third  point.  Objective  functions  for  a 
homogeneous  and  an  inhomogeneous  cost  surface  are  now 
considered . 

Figure  5  presents  the  elliptical  contours  (isodapanes) 


34 


for  an  objective  function  with  a  homogeneous  cost  surface 
C(x,y)=5.0.  This  cost  surface  is  plotted  in  Figure  1.  A 
point  anywhere  on  the  straight  line  between  the  two  end 
points  located  at  P ( - 2 , 2 )  and  P ( 2 , - 2 )  provides  an  optimum 
solution  for  the  two  link  path.  Substituting  the  values  in 
(17)  the  objective  function  becomes: 

F  ( x , y )  =5 . 0 1  x+2 )  2+  ( y-2  ) 2  +J(  x-2  )  2+  (y+2  ) 2 ) . (18) 

Using  the  property  of  distance  and  the  triangular 
inequality,  it  can  be  shown  that  this  objective  function  is 
convex.  The  objective  function  is  the  sum  of  two  distance 
functions  multiplied  by  a  constant.  The  distance  function  is 
convex  and  the  sum  of  two  distance  functions  multiplied  by  a 
constant  is  also  convex  indicating  a  convex  objective 
function.  The  particular  minimum  found  for  the  objective 
function  in  Figure  5  is  at  P ( x , y )=( 0 . 3696 0 . 3696 ) . 

Figure  6  demonstrates  the  contours  of  the  two  link 
objective  function  for  the  cost  surface  shown  in  Figure  3. 
The  cost  surface  is  given  by: 

C(  x,y)  =  10.  Oexp(- 0.2x2 -0.8y2)  +  1  >0 . (19) 

The  base  cost  is  incorporated  into  this  cost  surface  and  is 
given  the  value  of  1.0.  The  objective  function  now  becomes 

F(x,y)=Co(x,y)  (  ^/(x+2  )  2+(y-2  )  2  )  + 

Cf  ( x ,  y )  (  J(x- 2)*+(y+2)2  ) . (20) 

where  Co(x,y)  and  Cf(x,y)  are  the  average  costs  of  the  first 
and  last  link  evaluated  using  the  cost  surface  described  by 


35 


Figure  5  Contours  and  two  link  path  for  objective  function 
(18) 
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Fiaure  6  Contours  and  two  link  path  for  objective  function 
(20) 
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(19)  above.  The  fixed  end  points  are  located  at  P(-2,2)  and 
P ( 2 , -  2 )  .  A  minimum  value  of  the  objective  function  is  found 
at  P ( x , y ) = ( -4 . 1 1 29 , -2 . 0 1 54 ) .  The  Hessian  evaluated 
numerically  at  this  point  is: 

0.8174378  0.  1  114070 
0. 1 114070  2.2476149. 

The  eigenvalues  for  this  Hessian  are  0.8088115,  2.256241, 
indicating  that  the  matrix  is  positive  definite  and  that  the 
minimum  found  is  a  local  minimum.  Because  the  cost  surface 
is  symmetric  about  the  axes,  another  local  minimum  can  be 
found  in  the  other  pit  shown  by  the  contours.  It  is 
speculated  that  these  two  local  minima  belong  to  the  set  of 
global  minirrfa  for  this  objective  function.  For  large  |x|  or 
| y i  the  cost  surface  approaches  the  value  of  1.0.  In  this 
region  of  the  plane  this  particular  objective  function 
behaves  like  the  distance  function  with  increasing  value.  A 
minimum  in  this  region  is  unlikely. 

Clearly,  the  optimal  solutions  for  the  two  link 
objective  functions  with  both  homogeneous  and  inhomogeneous 
cost  surfaces  can  be  found  graphically.  Such  is  not  the  case 
for  paths  with  more  than  one  locational  variable.  The 
objective  function  surface  for  these  paths  extends  into  a 
space  of  more  than  the  three  dimensions  which  is  not  easily 
visualized.  The  contouring  of  such  surfaces  is  virtually 
impossible  and  the  numerical  solution  techniques  outlined 
above  are  required. 
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3.2.2  Variable  x  and  y  decision  variables 

An  obvious  choice  for  the  decision  variables  are  the  x 
and  y  coordinates  of  the  locational  variables  as  described 
in  the  two  link  method  above.  By  allowing  each  (x,y)  pair  to 
locate  freely  on  the  geographic  plane  one  can  argue  that 
some  location  of  such  pairs  should  satisfy  the  objective 
function  and  produce  a  minimum  cost  path.  The  objective 
function  using  these  decision  variables  is  given  by: 

n-l 

F(X,Y)=V(Ca(x.  ,  y.  )+BC)(  J(x.,,  -x.  )2+(y.  -y  )2  . (21) 

i  i  %  i+l  i  i+i  i 

i -0 

where  n  is  the  number  of  links  in  the  path,  x  and  y  are  the 

coordinates,  X=x,  , . . . ,x  ,  Y=y. , . . . ,y  ,  and  BC  is  the  base 

1  n-i  1  n-l 

cost.  Ca  equals  the  average  cost  of  link  i  and  is  given  by: 
m 

VC(x..  ,y..  )/m,  where  x..  =(x.  -x.  )j/m  and 

jr;  ij  »J  U  '+1  i 

y..  =  (y.  -y  )j/m.  For  a  n  link  path  2n  decision  variables 

ij  i+l  i 

are  required. 

The  x,y  decision  variable  method  is  illustrated  in 
Figure  7.  The  cost  surface  and  end  points  are  the  same  as 
for  Figure  6.  The  number  of  links  n  in  this  path  is  8  and 
the  cost  of  each  link  is  averaged  at  m=5.  The  number  of 
decision  variables  is  14.  The  solution  found  in  Figure  7  is 
a  local  minimum  because  the  eigenvalues  of  the  Hessian  are 
al 1  posi tive. 

3.2.3  The  variable  r  and  theta  method 

A  second  method  for  defining  the  decision  variables  of 
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F igure 


ClX.Y)  PARAMETERS: 
K  -  10.00 

Ar  0.20 

B-  0.80 

PATH«  1 . 

BC=  1.00 
C0ST=  12-35 
D I ST=  10.44 


^^-2.2) 


7  A  path  calculated  using  the  x,y  method 
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the  objective  function  is  proposed.  Let  the  path  be  defined 
by  a  fixed  number  of  equal  but  variable  length  links.  Let 
each  of  these  links  vary  their  orientation  on  the  x,y  plane. 
Then  the  decision  variables  are  the  angles  "theta"  each  link 
makes  with  the  x-axis,  and  the  length  "r"  of  each  and  every 
link.  Using  this  method  the  path  is  represented  by  the 
following  schema: 


Figure  8  A  five  link  path  with  equal  but  variable  link 
length 


The  decison  variables  for  this  representation  are  the  angles 
theta  tj  ,  .  .  .  ,  t4  and  the  link  length  r.  The  coordinates  x 
and  y  of  the  locational  variables  P(x,y)  are  given  by: 


xi  =\,+r  £  cos  ^ 

i=i 

i 

y  =y„  t  £  sin  t: 
'  i=i  1 


(22) 
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where  n  equals  the  number  of  links. 

In  order  to  incorporate  the  decision  variables  t,  and  r 
into  the  objective  function  special  attention  is  given  to 
the  relationship  between  the  last  link  and  r.  N  steps  of  r 
and  respecive  angles  t  defines  the  location  of  all  links  of 
the  path,  including  the  last.  If  the  path  is  to  reach  the 
end  point  with  last  link  of  length  r  exactly,  then  r  must  be 
set  to  the  length  of  the  last  link.  This  is  achieved  through 
the  use  of  a  quadratic  penalty  function  added  to  the 
objective  function.  The  minimum  of  the  squared  difference 
between  the  length  of  r  and  the  length  of  the  last  link 
occurs  when  the  two  lengths  are  equal,  and  the  contribution 
of  the  penalty  function  to  the  total  value  of  the  objective 
function  is  zero.  This  method  is  theoretically  sound,  but  it 
has  a  few  practical  problems.  As  the  optimization  proceeds, 
the  length  of  r  is  acted  upon  by  two  forces.  The  optimum 
region  of  the  path  requires  one  length  of  r  and  the  penalty 
function  requires  another.  When  the  path  is  located  near  the 
optimum,  small  changes  in  r  required  to  zero  the  penalty 
function  do  not  cause  large  changes  in  the  value  of  the 
objective  function  and  often  a  stationary  point  resulted. 
Fortunately,  the  solution  at  this  stationary  point  is  close 
to  a  local  minimum  because  the  gradient  is  near  zero. 

The  variable  link  length  method  discussed  above  is 
incorporated  into  the  following  objective  function: 
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n-1 

F  (  T  ,  r  )  =  ^^(  Ca  (  xj  ,  ^  )  +BC  )  r  + 

i  =0 

w(  yTxf  -x^)  2  +  (yf  -y  )  2  -  r)2  ,  . (23) 

where  n  is  the  number  of  links  in  the  path,  Ca  is  as  defined 
in  Section  3.2.2.,  r  is  the  link  length,  T=t,  ,...,t  are  the 
angles  each  point  makes  with  the  x-axis  (see  Figure  8),  BC 
is  base  cost  of  each  link  and  w  is  a  constant  weight.  A 
further  explanation  of  the  objective  function  (23)  is 
required.  The  last  term  is  actually  a  penalty  function  added 
to  the  cost  of  the  path.  For  a  large  enough  value  of  w,  the 
length  of  r  will  be  set  equal  to  the  length  of  the  last 
link.  In  order  that  the  penalty  function  does  not  conribute 
to  the  cost  part  of  the  objective  function,  the  last  link 
and  r  must  have  the  same  value.  Through  a  series  of 
experiments  a  weight  of  w=10.0  was  chosen.  This  weight 
brought  r  close  enough  to  the  length  of  the  last  link  such 
that  at  the  minimum,  the  penalty  function  contributed  very 
little  to  the  value  of  the  objective  function.  This  weight 
of  w  is  used  for  all  future  calculations  for  the  variable  r 
objective  function. 

Figure  9  illustrates  a  path  calculated  using  the 
variable  r  and  theta  method.  The  problem  posed  is  the  same 
as  that  for  Figure  7.  For  the  8  link  path,  the  variable  r 
and  theta  method  requires  8  decision  variables.  The  path 
demonstrates  the  equal  length  of  links.  The  path  in  Figure  9 
is  declared  a  local  minimum. 
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C(X.Y)  PARAMETERS: 
K=  10.00 

Ar  0.20 
Br  0.80 

PATH«  1  • 

BC=  1-00 
C0ST=  12.37 
D  I  ST=  10.45 


Figure  9  A  path  calculated  using  the  r  and  theta  method 
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3.2.4  A  comparison  of  the  two  methods 

The  variable  x,y,  and  r  and  theta  methods  can  best  be 
compared  by  two  numerical  experiments.  The  first  experiment 
is  designed  to  compare  the  location  of  the  path  calculated 
using  each  method.  Figure  10  illustrates  that  the  paths 
calculated  with  each  method  are  located  in  the  same  region 
of  the  plane.  The  only  noticable  difference  in  this  figure 
is  the  cost  contributed  to  each  path.  The  variable  x,y 
method  provides  a  path  with  0.02  less  cost.  This  difference 
in  cost  can  be  explained  by  the  path  shown  in  Figure  7, 
which  is  calculated  using  the  x,y  method.  In  the  highly 
curved  regions  of  the  path  the  link  length  is  shorter.  In 
the  more  straight  regions,  the  link  length  increases.  As  a 
result  the  x,y  method  produces  paths  which  are  smoother  in 
their  curved  extents.  This  smoothness  may  result  in  more 
cost  effective  paths.  The  second  experiment  illustrates  the 
performance  of  each  method  for  calculating  paths  with 
increasing  number  of  links.  Figure  11  clearly  indicates  that 
the  cost  of  the  path  for  each  method  converges  as  the  number 
of  links  increase.  At  10  links,  for  example,  the  cost 
becomes  virtually  identical.  However,  the  number  of 
iterations  required  to  calculate  the  path  for  increasing 
number  of  links  is  radically  different.  Figure  12  shows  that 
for  the  variable  x,y  method,  the  number  of  iterations 
increases  at  10  to  the  power  of  a  constant  times  the  number 
of  links.  The  variable  r  and  theta  method  converges  at  a 
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C(X.Y)  PARAMETERS: 

K=  10.00 
A=  0.20 
B=  0.80 

PATHc  1  . 

BC=  1-00 
C0ST=  12-31 
DISTr  10.43 


Figure  10  Paths  calculated  using  the  x,y  and  r,  theta 
methods 
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Figure  11  Cost  calculated  at  different  number  of  links 


METHODS  X.  Y  AND  R.  THETA  COMPARED 
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different  number  of  links 
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relatively  constant  number  of  iterations  regardless  of  how 
many  links  were  chosen  in  the  path. 

A  possible  explanation  for  the  difference  in  the  number 
of  iterations  required  for  convergence  in  each  of  the  two 
methods  is  found  in  the  relationship  of  the  decision 
variables  to  one  another.  In  the  variable  x,y  method,  the  x 
and  y  for  each  locational  variable  effects  only  the  adjacent 
links.  As  a  result  a  change  in  x  or  y  at  each  iteration 
effects  the  location  of  only  two  links.  Many  iterations  are 
required  to  transmit  the  change  down  all  the  links  of  the 
path.  In  the  variable  r  and  theta  method  a  change  in  the 
orientation  of  the  first  link  for  example,  effects  the 
location  of  the  last  link.  This  observation  is  supported  by 
the  relationship  defined  in  formulation  (22).  Because  the 
variable  r  and  theta  method  provides  virtually  the  same 
solution  for  paths  with  greater  than  9  links  and  for  fewer 
iterations,  it  appears  that  it  is  the  prefered  method  for 
calculating  minimum  cost  paths.  Other  numerical  experiments 
have  supported  this  observation. 

3.2.5  Numerical  integration  analysis 

The  objective  functions  (21)  and  (23)  are  in  part 
numerical  integrations.  Confidence  in  it's  ability  to 
approximate  a  continuous  integration  is  demonstrated  by 
Figure  13.  The  cost  surface  is  represented  by 
C ( x , y ) = 1 0 . Oexp ( -0 . 2x 2-0 . 8y2 ) .  Four  paths  between  the  same 
end  points  were  calculated  with  2,4,8  and  16  links.  The 
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C(X.Y)  PARAMETERS : 


Kr 

10.00 

A  = 

0.20 

B  = 

0.80 

PATH* 

1  . 

BC  = 

1  .00 

COST  = 

16.72 

0 1  ST  = 

9-59 

PATH* 

2. 

BC  = 

1  .00 

COST  = 

12.83 

D I  S  T  = 

10.S5 

PfrJH* 

3  • 

Bt  = 

1  .00 

COSf,= 

12.36 

DISTr 

10.43 

^4th* 

4  . 

BC  = 

1  .00 

COST  = 

12.27 

D I  ST  = 

10.41 

Figure  13  Paths  calculated  with  2,  4,  8,  and  16  links 


distance,  cost,  F  value,  number  of  iterations  and  CPU  time 
for  each  path  were  calculated  and  presented  in  Table  1. 
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Number 
of  1 i nks 

Path 

U 

D i stance 

Cost 
( area ) 

F 

value 

Number  of 
i terat ions 

CPU 

seconds 

1  2 

1 

9.594 

16.724 

16.841 

7 

0.054 

1  4 

2 

10.553 

12.826 

12.832 

19 

0.  140 

1  8 

3 

10.434 

12.359 

12.360 

24 

0.309 

1  16 

4 

10.411 

12.273 

12.273 

25 

0.491 

Table  1  Comparison  of  paths  with  varying  number  of  links 

The  paths  plotted  in  Figure  13  all  locate  within  close 
proximity  of  each  other.  The  paths  with  greater  than  two 
links  are  very  similar  in  both  cost  and  distance  as  can  be 
observed  in  Table  1.  The  difference  in  the  precision  of  the 
integration  between  paths  with  8  or  16  links  indicates  that 
the  choice  for  the  number  of  links  is  a  subjective  decision. 
Figure  11  illustrates  that  the  cost  approaches  an  asymptote 
as  the  number  of  links  increase.  A  path  calculated  with  10 
links  appears  to  fall  within  the  asymptotic  portion  of  the 
curve.  The  distance  and  cost  of  the  path  calculating  using 
Rankin's  (1979)  method  was  10.404  and  12.247  respectively. 
This  compares  favourably  with  the  16  link  path. 
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3.2.6  Other  solution  criteria 

A  solution  to  the  objective  functions  (21)  and  (23) 
must  satisfy  other  numeric  and  graphic  criteria  before  it 
can  be  accepted  as  a  reasonable  solution.  The  numeric 
criteria  requires  the  same  value  of  F  for  paths  with  the 
same  optimum  solution  and  different  starting  solutions;  and 
the  same  value  of  F  for  symmetrical  solutions.  Graphic 
criteria  is  found  in  the  location  of  the  path.  The  numeric 
criteria  may  indicate  an  optimal  solution,  but  the  solution 
may  not  make  any  graphic  sense.  All  solutions  must  be 
visually  examined.  Figure  14  shows  two  different  starting 
solutions  which  converge  to  the  same  optimum  solution.  The 
cost  surface  is  given  by  C ( x , y )  =  1 0 . Oexp ( - 0 . 2x 2 - 0 . 8y 2 ) .  The 
paths  plotted  are  path  1  with  optimum  solution  path  2  and 
path  3  with  optimum  path  4.  Path  2  and  path  4  are  identical 
as  indicated  by  the  values  of  F  and  path  distance. 

Sometimes  there  may  be  more  than  one  equivalent 
solution  to  a  particular  problem.  Intuitively,  if  both  the 
cost  function  and  the  beginning  and  end  points  are  symmetric 
about  an  axis,  there  should  be  two  symmetrical  solutions. 
Figure  15  demonstrates  two  virtually  identical  solutions  for 
a  cost  surface  given  by  C ( x , y ) = 1 0 . Oexp( -0 . 9x2-0 . 1y2 ) . 
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C1X.Y)  PARAMETERS'- 


K=  10.00 
A=  0.20 
B=  0.80 

PATH*  1  . 

BC=  1.00 
C0ST=  20.53 
D  I  ST=  11.24 


PATH" 


2  . 

1  .00 

16.74 

7.90 

3  . 

1  .00 

35.39 

4.02 

4  . 

1  .00 

16.74 

7.90 


Figure  14  Two  initial  solutions  -  one  optimum  solution 
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Figure  15  Two  symmetrical  solutions 
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3.3  Constrained  solution  paths 

Not  all  minimum  path  problems  can  be  expressed  by  the 
objective  functions  (21)  and  (23).  Some  minimum  path 
problems  require  that  the  solution  path  be  restricted  to 
certain  areas  of  the  geographic  plane  either  by  natural  or 
man  imposed  barriers.  Other  paths  may  be  subjected  to 
penalties  imposed  by  the  geographic  space  they  traverse.  In 
some  cases  such  problems  can  be  formulated  by  subjecting 
these  objective  functions  to  constraints.  These  constraints 
are  expressed  in  two  ways:  first,  through  the  use  of  linear 
constraints,  and  secondly,  through  the  use  of  penalty 
funct ions . 

Linear  constraints  formulated  by  equation  (15)  are  used 
for  finding  minimum  paths  constrained  by  linear  barriers  on 
the  x,y  plane.  For  each  locational  variable  there  exists  a 
region  beyond  which  the  constraints  are  violated.  Inside 
this  region  the  minimization  proceeds  as  if  no  constraints 
exist.  Figure  25  in  Chapter  4  demonstrates  a  constrained 
path.  For  efficient  optimization,  starting  solutions  within 
the  feasible  region  of  the  surface  are  preferred.  These  are 
Known  as  feasible  solutions.  Hauer  (1974)  expresses  linear 
inequality  constraints  as  a  linear  system  described  by 
G(X)=AX  +  B<0.  The  matrix  A  and  vector  B,  and  the  decision 
variables  X  representing  the  x,y  coordinates  of  the  path, 
can  force  the  path  to  be  restricted  to  certain  regions  of 
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the  plane. 

Constraining  the  path  by  adding  an  artificially  high 
cost  to  the  cost  surface  in  the  region  of  interest  may  be 
achieved  through  the  use  of  penalty  functions.  This 
technique  suggested  by  Werner  (1968)  is  concerned  with  using 
a  function  to  force  a  particular  locational  variable  to 
locate  outside  the  range  of  influence.  The  use  of  penalty 
functions  can  be  incorporated  into  any  cost  surface  and  no 
special  optimization  procedure  is  required. 


Chapter  4  Examples  of  Minimum  Cost  Paths 


The  location  of  minimum  cost  paths  is  as  conjectured  in 
previous  chapters  generally  not  achieved  through  intuitive 
insights.  By  presenting  examples  of  minimum  cost  paths  this 
chapter  demonstrates  two  points.  First,  that  the  problem  of 
where  to  locate  minimum  cost  paths  is  not  trivial  and 
second,  that  the  method  can  be  used  to  solve  different  types 
of  minimum  cost  path  problems. 

The  examples  presented  in  this  chapter  are  divided  into 
three  sections.  Section  4.1  presents  a  variety  of  problems 
exploring  some  general  relationships  between  the  cost 
surface,  the  end  points  and  the  minimum  cost  path.  The 
second  section  concentrates  on  finding  minimum  cost  paths 
for  cost  surfaces  which  are  derived  from  the  current 
literature.  By  providing  identical  examples,  the  methods 
used  in  this  thesis  are  supported.  Section  4.3  explores  the 
effect  of  simple  constraints  on  path  location.  The  examples 
which  use  constraints  further  support  the  use  of 
mathematical  optimization.  In  all  examples  presented  the 
number  of  links  was  kept  at  10  and  the  penalty  weight  w  of 
the  objective  function  was  set  at  10.0. 

4.1  Variations  on  one  cost  surface 

During  the  testing  phase  of  the  mathematical 
optimization  an  elliptical  negative  exponential  cost  surface 
provided  many  characteristics  of  a  generalized  cost  surface. 
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This  same  surface  is  used  in  this  section  to  demonstrate  the 
effect  of  cost  on  path  location.  The  surface  is  defined  by: 

C ( x , y ) =Kexp( -Ax2-By2 ) ,  . (24) 

where  K , A ,  and  B  are  positive  constants.  The  behaviour  of 
this  surface  is  ideal  for  presenting  problems  within  this 
section.  At  C(x,y)  =  K  at  x=y=0.  For  positive  A  and  B  with 
large  |x|  or  large  |y|;  C(x,y)  approaches  0.  The  cost 
function  C(x,y)  is  defined  everywhere  in  the  plane  implying 
that  a  path  can  locate  anywhere  without  creating 
discontinuity  problems. 

4.1.1  Cost  surface  orientation 

Figure  16  presents  four  paths,  the  end  points  of  which 
are  located  along  the  x-axis.  The  major  axis  of  the 
elliptical  cost  surface  is  located  on  the  y-axis.  The 
ability  of  a  minimum  cost  path  to  avoid  the  higher  cost 
areas  is  diminished  as  the  end  points  approach  the  center  of 
the  cost  surface.  To  avoid  the  high  cost  center,  the  path 
must  initially  traverse  through  the  high  cost  region.  In 
other  words  it  sometimes  is  cheaper  to  traverse  a  high  cost 
region  than  to  avoid  it.  The  cost  versus  distance  plot  in 
Figure  17  illustrates  the  cost  profile  of  the  path.  In 
regions  outside  the  high  cost  areas  of  the  cost  surface  the 
paths  are  virtually  straight  lines  due  to  the  nearly 
homogeneous  nature  of  the  surface. 
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C ( X  .  Y  )  PARAMETERS : 

K=  10-00 

A=  0.80 

B=  0.20 


* 


* 


* 


x 


* 


* 


Figure  16  Four  paths  traversing  the  cost  surface 
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Path 

U 

Base 

Cost 

Cost 

Distance 

Cost/Di  stance 

Sinuosi ty 
Ratio 

2 

6.0 

43.59 

4.00 

10.90 

1.00 

1 

6.0 

64.55 

9.28 

6.95 

1  .  16 

3 

6.0 

84.95 

13.23 

6.42 

1.10 

4 

6.0 

106.74 

17.11 

6.24 

1  .07 

Table  2  Cost  and  distance  analysis  for  Figure  16 


Table  2  provides  numerical  information  for  cost  and 
distance  variables  concerning  the  paths  in  Figure  16.  For 
each  minimum  cost  path  the  base  cost,  the  total  cost  (area), 
distance  and  the  two  ratios,  cost  divided  by  distance  and 
sinuosity,  are  given.  The  sinuosity  ratio  is  a  measure  of 
path  curvature  and  is  calclated  by  dividing  the  path 
distance  by  the  straight  line  distance  between  the  two  end 
points.  The  ratio  is  always  greater  than  or  equal  to  1.0. 
Table  2  illustrates  that  as  the  end  points  are  located 
further  away  from  the  high  cost  center,  the  cost/di  stance 
and  the  sinuosity  ratio  decrease.  Path  2,  which  cannot  avoid 
the  high  cost  of  the  center,  chooses  a  virtually  straight 
path  resulting  in  a  sinuosity  ratio  of  1.00.  This 
illustrates  the  fact  that  the  decision  of  whether  to  go 
through  or  to  avoid  the  high  cost  regions  of  the  surface  is 
dependent  on  both  the  starting  and  ending  points  of  the  path 
and  the  shape  of  the  cost  surface. 

The  purpose  of  the  next  example  is  to  illustrate  the 
effect  that  the  orientation  of  the  cost  surface  has  on  path 
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location.  The  major  axis  of  the  surface  is  now  centered  on 
the  x-axis.  The  end  points  remain  the  same  as  in  Figure  16. 


Path 

# 

Base 

Cost 

Cost 

Di stance 

Cost/Di  stance 

Si nuosi ty 
Ratio 

2 

6.0 

47.46 

4.00 

11.86 

1  .00 

3 

6.0 

36.95 

6.01 

6.  15 

1  .50 

1 

6.0 

42.11 

9.59 

4.39 

1  .20 

4 

6.0 

53.83 

13.07 

4.  12 

1  .09 

5 

6.0 

68.  14 

16.80 

4.06 

1  .05 

Table  3  Cost  and  distance  analysis  for  Figure  18 


The  effect  of  the  orientation  is  illustrated  by 
comparing  paths  in  Figures  16  and  18.  In  both  cases  the 
cost/distance  and  the  sinuosity  ratio  decrease  as  the  end 
points  locate  away  from  the  center.  However,  the 
relationship  between  the  end  points  and  the  high  cost 
regions  result  in  different  cost/di  stance  profiles  as 
demonstrated  in  Figures  17  and  19.  The  most  illustrative 
profile  is  that  for  path  2  (Figure  16)  shown  in  Figure  17 
and  path  3  (Figure  18)  shown  in  Figure  19.  Both  paths  have 
the  same  end  points.  Path  2  is  located  directly  through  the 
center.  No  effort  is  made  to  deviate  to  a  lower  cost  region. 
The  orientation  of  the  cost  surface  suggests  that  such  a 
deviation  would  be  of  little  benefit,  because  the  distance 
traveled  to  gain  this  region  is  too  large  and  too  costly. 
Path  3  in  Figure  18,  however,  does  deviate  to  a  lower  cost 
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Figure  1 


C(X.Y)  PARAMETERS: 


K  = 

10 

.00 

A  = 

o. 

20 

B  - 

o . 

80 

PATH  2 


* 


* 


8  Four  paths  avoiding  the  high  cost  regions 
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region  and  the  extra  distance  of  this  path  over  the  direct 
path  2  resulted  in  a  cost  saving  of  10.51.  The  orientation 
of  the  surface  allowed  the  higher  deviation  to  be 
prof i table. 

4.1.2  The  effect  of  cost  on  path  curvature 

The  decision  of  whether  to  go  around  a  high  cost  region 
of  the  cost  surface  or  through  it  must  be  made  based  on  the 
calculation  of  the  minimum  cost  path.  The  purpose  of  the 
paths  in  Figure  20  is  to  demonstrate  the  effect  of  the  cost 
on  path  curvature.  The  maximum  curvature  of  the  minimum  cost 
path  occurs  in  the  region  where  avoiding  high  cost  is  still 
possible.  The  cost/di  stance  profile  given  in  Figure  21 
indicates  the  relative  areas  occupied  by  each  path. 
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Figure  20  The  effect  of  cost  on  path  curvature 
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Figure  21  Profile  of  cost  versus  distance  from  one  end  point 
for  paths  in  Figure  20 
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Path 

U 

Base 

Cost 

Cost 

D i stance 

Cost/Di  stance 

Si nuosi ty 
Ratio 

2 

4.0 

16.71 

4.02 

4.16 

1  .01 

3 

4.0 

18.54 

4.  12 

4.50 

1.03 

1 

4.0 

22.50 

4.04 

5.11 

1.10 

4 

4.0 

28.86 

4.90 

5.83 

1  .22 

5 

4.0 

32.73 

5.10 

6.42 

1  .28 

6 

4.0 

34.68 

4.65 

7.46 

1  .  16 

7 

4.0 

35.44 

4.07 

8.71 

1  .02 

8 

4.0 

35.59 

4.00 

8.90 

1.00 

Table  4  Cost  and  distance  analysis  for  Figure  20 


The  curvature  of  the  paths  is  given  by  the  sinuosity 
ratio  in  Table  4.  Some  paths  have  a  low  sinuosity  ratio 
because  they  are  in  a  nearly  homogeneous  region  of  the  cost 
plane.  Others  have  a  lower  sinuosity  ratio  because  they 
cannot  deviate  to  the  low  cost  regions.  This  observation  can 
be  made  from  paths  2  and  8.  These  paths  approach  a  straight 
line  even  though  they  are  in  different  regions  of  the  cost 
surface.  The  effect  of  the  cost  surface  in  path  2  is 
negligible  as  it  contributes  only  0.63  to  the  total  cost. 
Path  8  cannot  avoid  the  high  cost  region  of  the  surface  and 
chooses  a  path  which  deviates  only  0.02  coordinate  distance 
at  x=0  from  a  straight  line.  Path  5,  the  longest  path,  has 
end  points  the  position  of  which  allows  for  the  avoidance  of 
the  high  cost  areas. 
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4.1.3  Varying  the  base  cost 

One  of  the  more  interesting  examples  which  can  be  used 
to  illustrate  the  relationship  between  cost  (area)  and 
distance  is  provided  by  Figure  22.  The  base  cost  is 
increased  progress i ve 1 y  until  it  dominates  the  cost  surface. 
The  example  verifies  the  obvious.  A  million  dollar  per  mile 
divided  highway  is  less  effected  by  variations  in  local 
costs  than  a  small  country  road  which  winds  around  every 
slough.  The  higher  the  base  cost  the  straighter  the  minimum 
cost  path,  ceter i s  par i bus . 


Path 

n 

Base 

Cost 

Cost 

Di stance 

Cost/Di  stance 

Sinuosi ty 
Rat  io 

2 

1  .0 

12.31 

10.42 

1  .  18 

1  .84 

3 

2.0 

o 

o 

CM 

CM 

9.09 

2.42 

1.61 

4 

4.0 

38.70 

7.77 

4.98 

1.37 

1 

8.0 

67.04 

6.56 

10.22 

1  .  16 

5 

16.0 

114.86 

5.69 

20.19 

1  .01 

Table  5  Cost  and  distance  analysis  for  Figure  22 


The  major  observation  for  this  example  is  shown  by  the 
sinuosity  ratio.  As  the  base  cost  increases,  the  sinuosity 
ratio  decreases  indicating  that  a  route  with  higher  base 
cost  is  less  influenced  by  the  cost  surface. 

4.2  Radially  symmetric  cost  surfaces 

The  use  of  radially  symmetric  surfaces  to  represent 
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C(X.Y)  PRRBMETERS: 


K=  10.00 

n=  0.20 

B=  0.80 


PATH# 
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D 1  ST  = 


1  . 
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67.04 

6.56 

2. 

1  .00 

12.32 

10.43 

3. 

2.00 

22.00 

9.09 

4. 

4.00 

38.70 
7  .77 

5. 

16.00 

114.86 

5.69 


Figure  22  The  increasing  base  cost  (BC)  example 
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Figure  23  Profile  of  cost  versus  distance  from  one  end  point 
for  paths  in  Figure  22 
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some  cost  variable  within  an  urban  setting  has  been  given 
recent  attention  by  Wardrop  (1969),  Angel  and  Hyman 
(1970,1972,1976),  Z i tron ( 1 974 ) ,  and  Puu ( 1 978a , 1 978b ) .  The 
cost  variable,  usually  measured  as  the  inverse  of  velocity, 
is  a  function  of  the  distance  r  from  the  city  center.  In 
order  to  provide  a  model  of  velocity  that  is  both 
mathematically  expedient  and  functional  within  the  urban 
setting,  the  authors  have  assumed  that  the  velocity  is 
distributed  symmetrically  from  the  center.  This  assumption 
may  not  be  valid  for  all  cities,  but  as  Angel  and  Hyman 
(1976)  point  out  it  is  useful  to  illustrate  a  methodology 
which  derives  minimum  cost  paths.  Presenting  radially 
symmetric  examples  in  this  section  has  one  purpose.  The 
examples  shown  are  well  Known  in  the  current  literature.  By 
comparing  the  paths  generated  by  the  mathematical 
optimization  to  those  found  in  the  literature,  the  merits  of 
both  methods  can  be  evaluated. 

4.2.1  Wardrop' s  velocity  surface 

One  of  the  radially  symmetric  surfaces  used  by  Wardrop 
(1969)  varies  the  velocity  as  the  distance  away  from  the 
center.  The  velocity  is  given  by  V(r)=Kr,  where  K  is  a 
constant  and  r  is  the  distance  from  the  city  center.  In  this 
thesis,  the  inverse  of  velocity  provides  a  suitable 
formulation  of  the  cost  (time)  surface,  which  becomes: 

C  ( x ,  y )  =  1  /Kr . (25) 

where  r  =  x/x2+ y2  .  In  order  to  duplicate  Wardrop' s  example,  a 
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partial  family  of  paths  radiating  from  the  point  P (  1  ,  0 )  is 
shown  in  Figure  24.  This  family  is  virtually  identical  to 
Wardrop' s .  Path  1  is  in  the  form  of  a  semi-circle  from  the 
point  P(-1,0)  to  P  (  1  ,  0  )  .  This  path  can  be  used  as  a  further 
validation  of  the  mathematical  optimization  method. 

4.2.2  Angel  and  Hyman's  velocity  surface 

The  velocity  surface  used  by  Angel  and  Hyman  (1976)  was 
obtained  from  an  empirical  study  of  travel  times  by  the 
SELNEC  Transportation  Study  (1968)  for  the  city  of 
Manchester  Great  Britain.  The  velocity,  expressed  as  a 
function  of  r,  the  distance  from  the  city  center,  is  given 
by: 

V  (  r  )  =  A-B*exp  (  -K*r  ) . (26) 

where  A,  B,  and  K  are  constants.  For  this  example  A=24.9, 

B  = 1 6 . 9  and  K  =  0.56.  The  travel  velocity  at  the  city  center  is 
8  miles  per  hour.  The  velocity  approaches  24.9  as  the 
distance  from  the  center  increases.  The  time  surface,  which 
is  used  to  represent  the  cost  surface  is  as  in  Wardrop' s 
example,  the  inverse  of  the  velocity. 

Angel  and  Hyman  (1976)  generate  a  family  of  minimum 
cost  paths  radiating  from  a  point,  and  in  this  thesis  the 
same  end  point  was  choosen.  The  family  of  minimum  cost  paths 
(Figure  25)  corresponds  to  the  paths  found  in  Angel  and 
Hyman's  work.  All  paths  except  the  one  from  the  center  curve 
away  from  the  high  cost  regions  of  the  cost  surface. 

Apart  from  presenting  their  family  of  minimum  cost 
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C(X.Y)  PARAMETERS : 

K=  1.00 

A=  0.00 


Figure  24  Minimum  time  paths  for  Wardrop' s  velocity  surface 
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C(X.Y)  PARAMETERS: 

K=  0.56 


Figure  25  Minimum  time  paths  for  Angel  and  Hyman's  velocity 
surface 
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paths,  neither  Wardrop  (1969)  and  Angel  and  Hyman  (1976) 
provide  much  discussion  on  the  path  locations.  This  thesis 
does  not  part  with  their  tradition.  The  paths  shown  in 
Figures  24  and  25  are  used  soley  to  demonstrate  that  the 
mathematical  optimization  method  can  duplicate  their  minimum 
cost  paths.  This  demonstration  is  significant,  for  the 
equivalence  of  the  minimum  cost  paths  illustrates  that 
different  methods  can  be  used  to  solve  the  problem.  In 
addition,  these  examples  provide  a  framework  for  discussion 
of  the  transformation  method  which  is  discussed  in  chapter 
5. 

4.3  Constrained  minimum  cost  paths 

Some  problems  require  that  minimum  cost  paths  be 
restricted  to  specific  regions  of  the  geographic  plane.  For 
this  type  of  problem,  the  use  of  constraints  to  locate  such 
minimum  cost  paths  seems  appropriate.  This  section  presents 
examples  of  minimum  cost  paths  constrained  by  linear  and 
non-linear  constraints.  Minimum  cost  paths  restricted  to 
regions  of  the  plane  by  straight  line  boundaries  are  derived 
by  imposing  linear  constraints  on  the  objective  function. 

The  barrier  and  corridor  examples  given  below  demonstrate 
linear  constraints.  Non-linear  constraints,  invoked  by 
penalty  functions,  may  be  used  for  excluding  paths  from 
regions  which  can  not  be  demarcated  by  straight  lines.  The 
minimum  distance  from  a  point  constraint  is  provided  as  an 
example . 
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4.3.1  The  barrier  problem 

The  barrier  problem  is  best  presented  by  the  following 
hypothetical  case.  Given  a  cost  surface  within  two  countries 
divided  by  an  international  boundary  and  the  requirement 
that  all  paths  originating  and  terminating  in  one  country 
must  remain  in  that  country,  then  the  minimum  cost  path  is 
constrained  by  the  border .  The  example  in  Figure  26  involve 
two  countries,  say  Canada  and  the  U.S.A.,  with  a  cost 
surface  traversing  the  border  given  by  y=-1.  The  path  1 
constrained  by  the  international  border  is  1.1  times  as 
expensive  as  the  minimum  cost  path  2.  Although  no  cost 
surface  is  provided  by  the  road  map  of  British  Columbia,  the 
highway  between  Princeton  and  CranbrooK  is  a  potential  real 
world  example  of  such  a  constraint  in  action.  Highway  #3  is 
forced  to  locate  to  the  south  of  mountain  ranges  and  lakes 
but  is  constrained  to  stay  in  Canada  by  the  international 
border . 

4.3.2  Minimum  cost  paths  contained  within  a  corridor 

The  location  of  transportation  facilities  within  a 

corridor  is  common  in  both  urban  and  rural  settings. 
Constraints  on  the  dimensions  of  the  corridor  may  restrict 
the  optimum  location  of  minimum  cost  paths.  The  example  in 
Figure  27  illustrates  a  corridor  defined  within  two  straight 
lines  at  y=+1  and  y=-1.  As  in  the  previous  example  linear 
constraints  may  be  applied  to  the  objective  function  such 
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Figure  26  Path  constrained  by  a  barrier 
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Figure  27  Two  paths  contained  within  a  corridor 
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that  the  optimum  path  will  locate  within  the  corridor. 

4.3.3  Fixed  distance  penalty  functions 

Some  types  of  constraints  cannot  be  expressed  linearly 
and  may  be  incorporated  into  the  objective  function  through 
the  addition  of  penalty  functions.  In  practice  these  penalty 
functions  are  extreme  additions  to  cost  within  certain 
regions  of  the  plane  and  a  path  traversing  these  regions 
would  be  subjected  to  virtually  an  infinite  cost.  The 
example  chosen  assumes  two  points  and  a  path  with  the 
restriction  that  the  path  must  never  be  located  closer  than 
a  fixed  distance  (radius)  from  either  point.  This  situation 
could  arise  if  a  pipeline  was  not  allowed  to  locate  within  a 
given  distance  of  an  Eskimo  village.  The  penalty  functions 
are  applied  in  the  following  manner.  The  cost  is  added  to 
the  objective  function  by  the  penalty: 

exp ( 25 . 0 ( r-d ) ) ,  . (  27  ) 

where  r  is  the  radius  and  d  is  the  distance  of  the  path  from 
a  particular  point  (village).  If  (r-d)  is  positive  the  added 
cost  becomes  large  and  if  (r-d)  is  negative  the  added  cost 
is  virtually  zero.  Hence  the  path  would  locate  outside  the 
region  defined  by  the  point  and  the  radius. 

Figure  28  contains  a  cost  surface  with  two  constraints 
imposed.  Paths  1  and  2  are  calculated  with  a  base  cost  (BC) 
of  4.0  and  Paths  3,  4,  and  5  are  calculated  with  a  base  cost 
of  8.  Only  Path  5  is  unconstrained.  The  location  of  each 
path  is  dependent  on  the  starting  solution.  If  the  starting 
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Figure  28  Paths  constrained  by  circular  regions 


81 


solution  is  located  between  the  constraints,  then  a  local 
minimum  was  found  in  the  same  region.  Similarly,  if  the 
starting  solution  is  located  on  one  side  of  the  constraint, 
a  local  minimum  is  found  on  the  same  side.  Before  the 
optimal  path  can  be  declared,  both  local  minima  must  be 
found  and  compared.  Path  1  located  between  the  constraints 
cost  45.67,  while  Path  2  located  to  the  top  of  one 
constraint  cost  38.70.  In  this  case  Path  2  is  a  better  local 
minimum  and  therefore  the  optimal  path.  For  Paths  3  and  4 
the  situation  is  in  reverse.  Path  3,  between  the  constraints 
cost  69.11  and  Path  4  to  the  top  of  the  constraint  cost 
69.57.  The  unconstrained  optimal  path  is  located  as  in  Path 
5  and  a  similar  symmetric  path  exists  between  Paths  3  and  4. 
This  unconstrained  path  indicates  that  the  constraints  do 
restrict  the  location  of  the  optimal  path. 

The  location  of  a  minimum  cost  path  for  this  example 
cannot  be  left  entirely  to  the  optimization.  As  shown  the 
optimization  routines  converge  on  a  variety  of  paths 
depending  on  the  location  of  the  starting  solution  and  all 
combinations  of  paths  through  the  constrained  regions  must 
be  tested.  Only  three  feasible  combinations  of  paths  exist 
in  this  example;  around  each  of  the  constraints  and  between 
them.  In  figure  28  the  placement  of  the  constraints  with 
respect  to  the  end  points  and  the  cost  surface  is  symmetric, 
therefore,  only  two  combinations  of  paths  required  testing. 
Other  problems  may  require  the  use  of  both  combinatorial  and 
optimization  techniques  for  a  solution.  These  types  of 
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problems  may  be  topics  for  future  research. 
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where  ds=^(dx)2  +(dy)2  +(dz)2  .  This  problem  also  requires  a 
calculus  of  variations  approach  and  the  solution  to  the 
appropriate  Euler  equations  for  most  surfaces  cannot  be 
found  analytically.  For  general  cases,  some  numerical 
technique  must  be  used  to  evaluate  the  integral  and  the 
optimization  techniques  employed  in  this  thesis  are 
appropr i ate . 

The  integral  (28)  may  be  formulated  as: 


minimize  F= 


(y.  -y.  )2  +(z.  -z.  ) 


i  + 


i+i 


. . . (29) 


where  n  equals  the  number  of  links  in  the  discrete  path 
representing  s.  The  objective  function  is  optimized  using 
the  variable  r  method  and  is  represented  by: 


n-i 


F(T,r)=  XI  V2  +<Zi+l  'Zi*2  + 


i=0 


w 


(  f-xn_,  >2  -My,  -yn_!  )2  -  r 


(30 


where  z  =z(x.  ,  y.  )  and  z.  =z(x.  ,  y.  ).  The 

j+1  1+1  1+1  I  •  • 

representation  of  T,r,x,y,  and  w  is  defined  in  Section 
3.2.3.  Because  the  multiplier  effect  of  the  cost  term  is 
absent  from  the  objective  function,  its  optimization  is  less 
difficult  than  calculating  the  minimum  cost  path. 
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5.1.2  Geodesic  examples 

Great  circles  calculated  between  two  points  on  a  sphere 
are  geodesics.  Because  of  their  important  classical  role  in 
cartography,  the  first  example  illustrates  a  geodesic  on  a 
hemisphere.  The  second  example  demonstrates  the  geodesic  on 
a  cone . 

The  great  circle  curve  from  Edmonton  to  London  is  shown 
in  Figure  29.  The  contours  of  the  surface  correspond  to  the 
latitude.  The  projection  of  both  the  latitudes  and  the 
geodesic  on  the  equatorial  plane  is  analagous  to  an 
orthographic  projection.  (It  was  not  derived  by  conventional 
means,  using  the  transformation  equations.)  Both  Edmonton 
and  London  are  connected  to  the  pole  via  meridians  which  are 
also  great  circles.  According  to  the  Atlas  of  Alberta,  the 
actual  distance  between  Edmonton  and  London  is  6,796 
Kilometers.  The  geodesic  distance  calculated  by  the 
optimization  routines  is  6,825  Kilometers.  This  29  Kilometer 
difference  requires  an  explanation. 

A  number  of  errors  resulted  in  the  formulation  of  the 
Edmonton  to  London  great  circle.  The  x  and  y  coordinates 
which  are  located  on  the  equatorial  plane  were  rounded  to 
the  nearest  hundred  Kilometers.  The  second  problem  concerned 
representing  the  radius  of  the  earth.  The  earth  is  a  geoid 
which  bulges  at  the  equator.  The  radius  calculated  as  the 
distance  from  the  equatorial  plane  to  the  pole  is  smaller 
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Figure  29  The  Edmonton  -  London  geodesic 
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than  that  on  the  equatorial  plane.  This  smaller  radius  would 
result  in  a  shorter  great  geoid  distance. 

The  geodesic  shown  in  Figure  30  is  calculated  from 


P(-1,-1)  to  P(+1,+1)  for  a  conic  surface  given  by: 


(31  ) 


Fortunately,  the  geodesics  for  this  cone  can  be  expressed 
analytically.  Appendix  I  contains  the  analytical  solution  to 
the  Euler  equations.  The  theoretical  distance  calculated 
analytically  in  this  example  is  3.901.  The  distance  of  3.899 
was  calculated  by  the  optimization  routines,  which 
demonstrates  their  preciseness.  Another  way  of  calculating 
the  distance  would  be  to  develop  the  cone  on  a  plane  and 
measure  the  straight  line  distance  (geodesic)  between  the 
two  points.  The  eigenvalues  of  the  Hessian  for  the  objective 
function  are  all  positive,  indicating  that  this  geodesic  is 
a  local  minimum. 

5.1.3  The  geodesic  and  minimum  cost  path  compared 

During  the  preliminary  stages  of  this  research  the 
minimum  cost  path  was  often  compared  to  the  geodesic.  Two 
reasons  invalidate  this  comparison.  First,  the  geodesic  is 
measured  in  terms  of  distance,  the  minimum  cost  path  in 
terms  of  area.  Secondly,  the  geodesic  is  located  on  a 
surface  whereas  the  minimum  cost  path  is  located  on  the  x,y 
plane.  A  demonstration  of  the  location  of  the  paths  relative 
to  the  surface  contours  can  be  made  by  projecting  the 
geodesic  onto  the  x,y  plane.  Figure  31  presents  the  minimum 
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Figure  30  Geodesic  on  a  cone 
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cost  path  (path  1)  and  the  geodesic  (path  2).  For  this 
particular  example  the  two  are  projected  differently  onto 
the  x,y  plane.  In  addition,  the  geodesic  is  about  three 
times  as  long  as  the  minimum  cost  path.  For  radially 
symmetric  surfaces  with  one  end  point  at  the  center  the 
projection  of  the  geodesic  onto  the  x,y  plane  corresponds  to 
the  minimum  cost  path  as  can  be  observed  in  Figure  32. 
Generally,  however,  the  projected  geodesic  does  not 
correspond  with  the  minimum  cost  path. 

5.1.4  Geodesics  as  minimum  cost  paths 

Some  types  of  minimum  path  problems  can  be  formulated 
to  minimize  distance  rather  than  cost.  In  order  to  give 
geographical  meaning  to  the  geodesic  within  this  context  the 
x,  y,  and  z  coordinates  must  all  be  of  the  same  physical 
dimension  and  scale.  This  condition  restricts  the  use  of  the 
distance  measure  as  a  substitute  for  cost  to  applications 
for  which  a  3-dimensional  surface  is  defined;  such  as 
calculating  the  great  circle  distance  on  the  globe  as 
demonstrated  above.  This  class  of  problem  is  significant 
from  a  purely  theoretical  point  of  view,  however,  and 
warrants  some  exploration. 

A  theoretical  example  of  this  problem  is  presented  in 
the  following  scenario.  Assume  a  sinuous  landscape  in  either 
the  x  or  the  y  direction.  The  problem  is  to  determine  the 
location  of  a  minimum  length  path  between  two  points  on  the 
surface.  Examples  of  such  problems  can  be  found  in  the 
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Figure  31  Geodesic  and  the  minimum  cost  path  compared 
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Figure  32  Geodesic  and  the  minimum  cost  path  compared 
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process  of  locating  new  pipelines  or  utility  corridors. 

Figure  33  demonstrates  three  geodesics  for  a  surface 
given  by: 

Z(x,y)=4sin(0.  1x)cos(y)+4 . (32) 

The  geodesic  represented  by  path  1  is  calculated  from  the 
point  P(x,y,z)=(-6,6, 1 .83)  to  (6, -4, 2. 52).  Path  2  is 
calculated  from  (-8, -5, 3. 19)  to  (8,5,4.81)  and  path  3  is 
calculated  from  (-8, -6, 1.24)  to  (-8,6,1.24).  Each  of  these 
geodesics  seems  to  locate  in  regions  of  the  surface  where 
the  change  in  the  z  coordinate  is  small.  This  is  especially 
noticable  in  path  3.  A  path  generated  above  a  straight  line 
on  the  x,y  plane  generates  a  distance  of  26.16,  whereas  the 
length  of  the  geodesic  is  22.10.  A  result  such  as  this 
illustrates  that  it  is  sometimes  shorter  to  walk  around  a 
hill,  then  over  i t . 

5.2  Transformat  ions ,  geodesics  and  minimum  cost  paths 

The  previous  section  demonstrated  the  geodesic  as  an 
alternative  to  the  minimum  cost  path.  For  some  problems, 
such  as  finding  great  circles  on  the  globe,  the  geodesic  may 
be  the  only  feasible  form  of  the  minimum  cost  path.  The 
natural  extension  to  this  observation  is  to  find  the 
conditions  for  which  the  geodesic  can  be  substituted  for  the 
minimum  cost  path.  The  observations  in  the  previous  section 
imply  that  the  geodesic  which  represents  length  must  also 
represent  cost.  The  surface,  on  which  the  geodesic  is 
calculated,  must  therefore  equate  cost  and  distance.  The 
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Figure  33  Three  geodesics  on  a  sinuous  surface 
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type  of  cost  surface  considered  by  this  thesis  does  not  meet 
this  condition.  The  z  values  of  this  surface  indicate  cost 
at  some  coordinate  location  x  and  y.  If  this  cost  surface  is 
to  be  used  to  derive  the  cost  geodesic,  all  x,y,  and  z 
coordinates  must  depict  locations  in  the  cost  space.  This 
implies  that  the  cost  surface  must  be  transformed  to  a 
three-dimensional  cost  space  if  it  is  to  be  used  for 
calculating  geodesics. 

The  purpose  of  this  section  is  to  summarize  and  discuss 
the  results  of  the  tr ansformat i on  research  of  Warntz  (1968), 
Angel  and  Hyman  (1972a,  1972b,  1976),  and  Puu  (1977,  1978a, 
1978b).  Their  research  addresses  the  problem  of  transforming 
the  cost  surface  into  another  surface  where  upon  the 
geodesic  corresponds  to  the  minimum  cost  path. 

5.2.1  Puu' s  contribution 

The  theme  of  the  transformat  ion  method  is  presented  by 
Puu  ( 1978b) : 

The  concept  of  a  cost  surface 

The  optimal  paths  in  general  turn  out  to  be 
curved.  Owing  to  this  the  natural  question  has  been 
posed  whether  it  is  possible  to  map  the  original 
region,  R,  onto  some  other  region  S  in  such  a  way 
that  the  curved  paths  are  mapped  onto  straight  lines 
in  some  sense  and  so  that  distance  in  the  image 
region  equals  cost. 

Wardrop  (1969)  has  investigated  the 
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possibilities  of  mapping  R  onto  some  other  plane 
region  S  by  conformal  mappings,  that  is,  by  complex 
analytic  functions  with  nonvanishing  derivatives. 
Angel  and  Hyman  (1970,  1972,  1976)  on  the  other  hand, 
guided  by  conjecture  by  Warntz  (1967),  explore  the 
possibilities  of  mapping  R  onto  a  curved  surface  S, 
embedded  in  three-dimensional  euclidean  space,  where 
the  images  of  the  optimal  paths  are  'straight'  in  the 
sense  of  being  geodesic.  Their  discussion  too  is 
confined  to  conformal  maps,  characterized  by  the 
facts  that  angles  are  preserved  and  that 
magnification  is  independent  of  direction. 

The  reasons  why  discussions  are  confined  to 
conformal  maps  are  never  stated,  but  I  am  going  to 
demonstrate  that,  in  fact,  conformal  mapping  is  the 
only  type  that  works  with  isotropic  transportation 
problems.  Angel  and  Hyman,  moreover,  only  discuss  the 
special  case  of  surfaces  of  revolution.  They  suggest 
that  the  method  has  a  much  greater  degree  of 
generality  than  Wardrop's,  but  the  exact  field  of 
application  never  becomes  clear.  In  fact  the  general 
assumptions  they  use  turn  out  to  be  too  restrictive. 

Puu  then  continues  with  theorems  and  conclusions  to 
demonstrate  three  concepts  concerning  the  transformation 
me  t  hod . 

Puu  concludes  that  "the  cost  for  any  isotropic 
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transportation  problem  is  a  conformal  transformation  of  the 
plane".  An  isotropic  transportation  problem  is  a  model  in 
which  the  cost  of  tr anspor tat i on  or  movement  is  defined  at  a 
point.  The  cost  surfaces  explored  by  Warntz,  Wardrop,  Angel 
and  Hyman,  and  in  this  thesis  are  all  isotropic 
transportation  models.  Puu' s  transformed  cost  surface  is  the 
conformal  map  of  these  isotropic  models.  Next  Puu  develops 
the  necessary  conditions  which  must  be  met  by  the  cost 
surfaces  such  that  the  geodesic  on  the  transformed  cost 
surface  corresponds  to  the  minimum  cost  path.  He  verifies 
that  the  surfaces  used  by  Angel  and  Hyman  meet  the  necessary 
conditions  for  a  transformation  and  duplicates  their 
examples . 

Puu' s  final  conclusion  states  that  it  is  "always 
possible  to  find  a  (local)  cost  surface  such  that  distances 
equal  cost  in  an  isotropic  transpor tat  ion  problem".  He 
observes  that  from  these  local  effects  nothing  can  be 
concluded  about  the  global  nature  of  the  cost  surfaces  which 
generally  can  have  many  folds  and  complex  curvatures. 

5.2.2  The  three  steps  in  the  transformation  process 

The  actual  calculation  of  minimum  cost  paths  using  the 
transformation  methodology  is  a  three  step  process.  The 
first  step  requires  the  formulation  of  the  transformation 
which  maps  the  original  cost  surface  into  one  in  which 
distance  equals  cost.  The  second  step  requires  the 
calculation  of  the  geodesic  on  the  transformed  surface.  The 
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third  step  uses  the  reverse  transformat  ion  to  map  the 
geodesic  into  the  minimum  cost  path  on  the  original 
geographic  plane.  These  three  steps  warrant  some  discussion 
because  using  this  process  Angel  and  Hyman,  and  Puu  find 
minimum  cost  paths  for  only  two  cost  surfaces. 

The  conformal  transformations  for  a  specific  cost 
surface  where  cost  is  expressed  as  a  function  of  x  and  y  is 
in  general  not  easily  achieved.  Even  if  the  surface  has  this 
isotropic  characteristic,  the  transformation  requires  the 
application  of  Gaussian  differential  geometry.  According  to 
Tobler  (1961)  this  is  a  very  difficult  topic  and  best 
avoided.  These  difficulties  are  evident  in  the  work  of  Angel 
and  Hyman  (1972  b) ,  and  Puu  (1978  b).  Both  authors  have 
transformed  the  same  surfaces.  The  radially  symmetric 
velocity  surface  of  Wardrop  (1969),  which  is  demonstrated  in 
Section  4.2.1,  becomes  a  cylinder  when  transformed.  The 
velocity  surface  V(r)=Ar2  +  B,  where  A  and  B  are  constants, 
and  r  is  the  distance  from  the  city  center  is  transformed 
into  the  sphere.  Although  Puu  (1978b)  provides  the  necessary 
conditions  which  must  be  met  by  the  cost  surface  before  it 
can  be  transformed  from  the  region  R  to  the  region  S,  his 
examples  stop  at  the  radially  symmetric  cost  surfaces  of 
Angel  and  Hyman.  He  shows  that  these  surfaces  also  meet  the 
necessary  conditions  for  the  transformation,  but  indicates 
that  the  transformat  ions  themselves  are  difficult  to 
achieve . 

The  second  step  in  finding  the  minimum  cost  path 
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requires  that  the  geodesic  on  the  transformed  surface  be 
found.  The  methods  illustrated  in  Section  5.1  of  this  thesis 
can  be  used  once  an  analytic  form  of  the  surface  is  given. 
Angel  and  Hyman,  and  Puu  provide  analytic  surfaces  for  which 
the  geodesic  are  well  Known.  The  cylinder  produces  helices 
and  the  sphere  produces  great  circles.  Generally  however, 
the  geodesics  are  not  easily  found  by  intuitive  or  analytic 
methods.  Angel  and  Hyman  (1970)  use  Huygens  (1912)  method 
for  calculating  geodesics.  Isochrones,  lines  of  equal  value 
from  a  point,  are  constructed  using  the  original  cost 
surface.  The  family  of  geodesics  from  this  point  is  normal 
to  each  isochrone.  Angle  and  Hyman  (1972)  state  that 
Huygen' s  construction  requires  considerable  adaptation 
before  it  can  be  a  practical  method  for  computing  minimum 
cost  paths  for  general  surfaces. 

The  third  step  concerns  mapping  the  geodesic  back  into 
the  original  geographic  plane  from  the  cost  space  such  that 
a  minimum  cost  path  can  be  plotted.  This  mapping  depends  on 
the  initial  transformation:  if  it  exists,  then  an  inverse 
transformation  should  be  found  which  can  map  the  geodesic 
into  the  minimum  cost  path.  No  attempt  has  yet  been  made  to 
use  this  inverse  mapping  method. 

The  difficulties  in  applying  the  three  step 
transformat  ion  method  is  illustrated  by  Angel  and  Hyman's 
(1972  b)  velocity  surface.  Although  their  surface  is 
radially  symmetric  the  transformation  could  not  be  achieved 
analytically  and  required  numerical  methods.  Once  the 
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transformed  surface  was  derived,  the  geodesics  could  not  be 
found  and  the  transformat i on  was  used  only  to  illustrate  the 
form  of  the  time  surface.  Angel  and  Hyman  used  Huygen' s 
construction  to  plot  the  isochrones  for  this  velocity 
surface.  Minimum  cost  paths  as  shown  in  Figure  25  were  drawn 
perpendicular  to  these  isochrones. 

5.2.3  Results  of  the  transformation  methods 

One  result  of  Angel  and  Hyman's  (1972  b)  transformation 
process  is  the  proof  that  the  minimum  cost  path  of  the 
original  surface  corresponds  to  the  geodesic  on  the 
transformed  surface.  This  mathematical  proof  clearly  links 
the  minimum  cost  path  to  that  of  the  geodesic.  There  are 
theoretical  advantages  in  providing  this  link.  The  intuitive 
conjecture  of  Warntz  (1965)  was  shown  to  be  true.  "A 
transformation  exists  which  will  map  a  realistic 
distribution  of  transpor tat  ion  costs  on  the  Euclidean  plane 
into  a  curved  surface  with  uniform  transport  facility" 

(Angel  and  Hyman  (1972  b).  The  validation  of  this  conjecture 
provided  new  impetus  to  the  tr ansformat i on  methods.  Warntz' s 
example  of  minimum  land  acquisition  routes  in  the  United 
States  is  used  as  reference  in  much  geographic 
transportation  and  path  research.  According  to  Harvey  (1969) 
this  research  is  important  because  the  geographic  concept  of 
distance  can  no  longer  be  viewed  as  the  two  dimensional 
Euclidean  distance,  but  must  be  extended  to  more  dimensions 
to  take  into  account  the  activities  which  may  effect  this 
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distance.  The  transformation  methods  contribute  to  this 
exploration  of  the  distance  theory  in  geography. 

There  are  however  major  difficulties  in  pursuing  the 
transformation  methods  for  deriving  minimum  cost  paths.  The 
difficulty  lies  in  the  mathematical  conditions  which  must  be 
met  by  the  cost  surface.  Although  Puu  (1978  b)  provides  the 
necessary  conditions,  he  provides  no  example  beyond  those  of 
Wardrop  (1969),  and  Angel  and  Hyman  (1972  b).  Also,  no 
example  of  transforming  a  non-symmetr i c  cost  surface  is 
found.  Puu  indicates  that  transformation  can  always  take 
place  locally  for  an  isotropic  cost  surface,  but  he  suggests 
that  providing  the  actual  transformations  and  geodesic  is 
very  difficult.  Thus,  if  non-radially  symmetric  cost  surface 
such  as  the  ones  used  in  Chapter  4  are  to  be  transformed, 
major  mathematical  research  must  be  attempted. 


Chapter  6  Conclusion 


This  thesis  solves  two  problems  in  the  calculus  of 
variations;  the  minimum  cost  path  problem  and  the  problem  of 
finding  the  geodesic  on  a  given  surface.  For  each  of  these 
problems  an  objective  function  is  formulated  and  analyzed, 
and  a  number  of  case  studies  are  performed.  This  chapter 
summarizes  the  results  of  using  mathematical  optimization 
techniques  to  solve  these  problems. 

6.1  The  methodology 

Perhaps  the  most  significant  contribution  of  this 
thesis  to  minimum  path  research  is  found  in  the  methodology 
of  deriving  the  minimum  cost  path.  Although  the 
representat ion  of  the  cost  integral  by  an  objective  function 
is  not  new,  the  development  of  the  r  and  theta  method  and 
the  use  of  constraints  have  not  yet  been  explored  in  the 
literature.  This  section  discusses  the  merits  and 
disadvantages  of  using  the  mathematical  optimization  method 
to  solve  the  Euler-Lagrange  equations. 

6.1.1  The  method  works 

Representing  the  cost  distance  and  the  distance 
integral  by  an  objective  function  which  is  then  minimized  is 
a  technique  that  works.  Two  analyses  referred  to  in  this 
thesis  support  this  conclusion.  The  trajectory  method  of 
Rankin  (1979)  calculated  a  similar  minimum  path  to  that 
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calculated  with  the  variable  r  and  theta  method.  This 
comparison  is  presented  in  section  3.2.5.  The  second 
analysis  is  shown  in  section  5.1.2.  The  length  of  the 
geodesic  is  calculated  using  both  analytic  and  optimization 
methods.  The  results  differed  by  a  distance  of  0.002.  Thus 
two  independant  methods  were  used  to  verify  the 
representation  of  the  integrals  (3)  and  (28)  by  the 
objective  function  (23)  and  (30)  respectively. 

Confidence  in  the  variable  r  and  theta  objective 
function  provides  many  benefits  for  minimum  cost  path 
research.  Many  different  path  problems  can  be  posed  and 
solved  as  demonstrated  in  Chapter  4.  The  method  is  flexible 
and  can  deal  with  the  interaction  of  the  end  points,  cost 
surface  and  base  cost.  The  varying  base  cost  problem  shown 
in  section  4.1.3  is  an  example  of  this  interaction.  The 
current  path  literature  has  shown  families  of  minimum  paths 
for  radially  symmetric  surfaces,  but  has  not  presented 
problems  in  which  the  base  cost  interacts  with  the  cost 
surface.  In  addition,  the  method  verifies  the  obvious 
minimum  cost  paths,  such  as  the  straight  line  of  the 
homogeneous  cost  surface,  and  the  geodesic  lines  of 
longitude  from  Edmonton  to  the  north  pole  and  from  London  to 
the  pole. 

Minimizing  a  non-linear  objective  function  with  many 
variables  also  presents  some  problems  regarding  the 
acceptance  of  a  minimum  solution.  Although  there  are  strong 
indications  that  the  objective  functions  are  convex,  the 
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minimization  of  some  specific  problems  resulted  in  solutions 
which  could  not  be  accepted  as  optimum.  For  example,  the 
cost  surface  represented  by  C ( x , y ) =Kexp ( - Ax 2 -By 2 )  becomes 
virtually  zero  for  large  |x|  or  |y|.  The  derivatives  of 
C(x,y)  with  respect  to  x  or  y  in  this  region  are  also  near 
zero.  If  no  base  cost  is  added  to  this  function  a  path  has 
virtually  zero  cost  and  can  locate  anywhere  in  the  plane.  A 
zero  gradient  is  one  of  the  necessary  conditions  for 
accepting  a  minimum  cost  path.  Hauer's  (1974)  algorithm 
indicated  a  zero  gradient  for  solutions  with  large  |x|  or 
| y |  but  simple  intuitive  reasoning  indicated  that  such 
solutions  are  meaningless  in  reality.  The  solution  methods 
could  not  differentiate  between  a  zero  gradient  and  a  zero 
derivative  of  the  cost  surface.  The  Hessian,  however,  for 
this  case  is  not  positive  definite  and  no  local  minimum  is 
claimed.  The  addition  of  a  base-cost  resulted  in  more 
reasonable  minimum  paths. 

6.1.2  The  constraints 

The  derivation  of  the  minimum  path,  which  is  subjected 
to  some  type  of  locational  constraint  is  sometimes  very 
difficult.  Current  literature  avoids  the  discussion  of  this 
topic  entirely  and  the  examples  presented  in  this  thesis 
were  of  a  very  simple  exploratory  nature.  The  reason  for  the 
difficulty  lies  more  in  the  practical  than  in  the 
theoretical  nature  of  constraints.  It  is  easy  to  draw  a  line 
over  which  the  minimum  path  cannot  cross,  but  it  is  more 
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difficult  to  formulate  this  constraint  into  the  solution 
procedures  of  the  optimization.  These  difficulties  are 
overcome  for  two  types  of  constraints.  The  simple  linear 
inequality  constraint  and  the  representat ion  of  a 
constrained  region  by  penalty  functions.  The  examples 
presented  in  section  4.3  illustrate  initial  success  at  the 
use  of  constraints. 

The  locational  constraints  expressed  as  linear 
inequality  constraints  presented  the  fewest  problems  to 
Hauer's  (1974)  algorithms.  Provided  that  the  initial 
starting  solutions  were  within  the  feasible  or  unconstrained 
region,  the  solution  was  found  easily.  Difficulty,  however, 
was  encountered  in  applying  these  linear  constraints  to  the 
variable  r  and  theta  objective  function.  The  decision 
variables  for  this  method  are  not  easily  subjected  to 
constraints  imposed  on  the  location  of  the  path.  Instead  the 
variable  x,y  objective  function  formulation  is  used 

Some  locational  constraints  could  not  be  expressed  by 
linear  constraints  and  penalty  function  were  applied. 

Section  4.3.3  presented  an  example  in  which  two  circular 
regions  of  the  plane  were  subjected  to  a  cost  penalty.  The 
minimum  paths  were  located  outside  this  region  indicating 
that  the  penalty  function  worked.  There  is  however  a 
topological  problem  related  to  this  problem.  The  minimum 
path  is  found  by  an  iterative  procedure  which  depends  upon 
the  starting  solution.  A  starting  solution  between  the  two 
penalized  regions  resulted  in  a  minimum  path  also  located  in 
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the  same  region.  Similarly,  if  the  starting  solution  is 
located  on  one  or  the  other  side  of  the  penalized  region, 
then  the  minimum  path  is  located  within  the  same  region. 
Conceivably,  if  the  number  of  such  penalized  regions  is 
large,  the  number  of  solution  paths  which  must  be  tested  is 
dependent  on  the  number  of  possible  feasible  paths  through 
these  regions.  This  is  a  combinatorial  problem  and  each 
feasible  path  must  be  derived  from  the  optimization. 

However,  given  a  particular  starting  solution  and  the  use  of 
a  penalty  function,  the  minimum  path  reflected  the  topology 
of  the  starting  solution.  The  acceptance  of  this  solution  as 
the  global  optimum  one  can  not  be  made.  Only  if  all  possible 
solution  topologies  are  tested  could  such  a  decision  be 
reached . 

6.2  The  transformation  method 

The  limited  analysis  and  discussion  of  the 
transformation  methods  given  in  this  thesis  is  justified  for 
two  reasons.  One,  the  methods  are  difficult  in  terms  of  the 
explanatory  mathematics  required  and  further  development  of 
these  methods  is  out  of  the  scope  of  this  thesis.  Two,  the 
ability  of  the  transformat  ion  method  to  solve  the  variety  of 
problems  posed  in  Chapter  4  of  this  thesis  has  not  been 
shown,  indicating  that  it  is  still  in  a  development  stage. 
The  mention  of  these  methods  in  Chapter  5  is  however 
relevent,  because  the  calculation  of  the  geodesic  is  a  topic 
that  is  discussed  in  the  current  transformat  ion  literature. 
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6.3  Future  research 

6.3.1  The  cost  surface 

This  thesis  represents  cost  by  a  continuous  function  of 
x  and  y.  Angel  and  Hyman  (1972  b)  provide  a  rough  fit  for 
their  negative  exponential  surface  and  refer  to  the  work  of 
Clark  (1951)  for  further  support.  It  is  recognized  that  the 
real  cost  space  is  in  general  very  complex  and  any  surface 
representing  this  space  may  contain  many  i r regu 1 ar i t i es .  In 
order  to  satisfy  the  continuity  condition  required  by  the 
optimization  procedures  only  limited  options  are  open  for 
representing  realistic  cost  surfaces. 

A  piecewise  continuous  cost  surface  representing  local 
regions  may  possibly  be  worked  into  the  solution  procedures. 
A  redesign  of  the  decision  variables  in  the  objective 
function  may  be  in  order  for  this  approach.  Finally  a  three 
dimensional  representation  of  cost  represented  by  C(x,y,z) 
is  explored  by  Werner  and  Boukidis  (1963)  and  Steenbrink 
(1974).  The  type  of  cost  surface  is  very  realistic  for 
mountainous  environments.  Again  the  decision  variables  of 
the  objective  function  must  be  modified  or  expanded  to 
include  the  third  dimension. 

6.3.2  The  use  of  constraints 

Two  immediate  projects  concerning  constraints  can 
result  from  this  thesis.  The  first  project  could  consist  of 
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expanding  the  use  of  linear  constraints  beyond  those 
presented  in  Sections  4.3.1  and  4.3.2.  The  second  project 
could  investigate  the  general  use  of  penalty  functions  to 
represent  a  wider  range  of  non-linear  constraints  than  is 
given  in  Section  4.3.3.  The  topological  problem  discussed 
earlier  may  be  part  of  this  project. 

6.3.3  Other  path  problems 

Chapter  4  and  5  demonstrated  the  feasibility  of  using 
objective  functions  to  find  minimum  cost  paths  and 
geodesics.  By  modifying  these  functions  a  variety  of  other 
path  related  problems  may  be  solved.  Typically,  these 
problems  may  deal  with  the  minimization  or  maximization  of 
some  interaction  between  a  path  and  a  surface.  The  attention 
given  to  such  path  problems  in  current  geographic  literature 
indicates  that  both  the  problems  and  their  solution 
methodologies  require  further  attention.  The  use  of  the 
mathematical  optimization  methodology  to  solve  these 
problems  is  as  shown  in  this  thesis  both  realistic  and 
feas i b le . 
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Find  the  Geodesic  on  a  cone  given  by: 


z(  x  ,y)  =  1  .  5y  x 2+y 2  . 

In  spherical  coordinates: 
x=r  sin(9)  cos(0) 
y=r  si n ( 9 )  sin  0 
z=r  cos(0) 

For  a  cone  9  =  a  constant 

Let  sin(0)=LY  and  cos(9)=/3  ,  then: 

x=  C^r  cos(0)  and  dx=Qf(cos(0)  dr  -  r  sin(0)  d0 ) 

y=C^r  sin(0f)  and  dy=Q'(cos(0)  dr  -  r  sin(0)  d0) 

z=  /3  r  and  dz=  /3  dr 

If  the  distance  of  a  geodesic  is  s;  then 
ds2  =  dx2  +  dy2  +  dz2 , 

which  implies  ds=(Q'2  r  2+  (dr/d0)  2 )  "  dO . 


In  order  to  minimize  the  integral 


a 


from  point  a  to  b,  the  following  Euler  equation  must  be 
solved : 


(33) 


(34) 


(35) 


dr/d0=O'r(  (r2/ro  2)-1  )  , 


(36) 


1  16 


where  r0  is  a  constant. 

Set  r/rQ  =sec(u)  then, 

dr/d0  =  ro  sec ( u ) tan ( u ) ( du/d0 ) 

and  CYr((r2/r0  2)-1  =  £Yr  sec  (  u  )  tan  (  u  )  , 

which  implies  r  =  rQ  sec{  CX  ( 0-0o  ) }  is  a  solution  to  the  Euler 
equation  (36).  If  r  is  substituted  into  (34)  and  integrated 
the  following  results: 

ds=(Qf2  r2  +  (dr/d(0)  2 )  60, 

dr  / d0=Q'  r  cosec 0-0o  )  }cot  {Ctf 0-0o  ) }  , 

and  ds=Q2  rQ  2  { sec2  ( u  ) +sec2  ( u )  tan 2  ( u  ) }  60, 

-(yr o  sec2  ( u  )d0; 


which  implies  that  s=£V  r0  /  sec2(u)d0, 


=  r0  tan{C^(  0-0o  ) } 


(37) 


In  order  to  calculate  the  length  of  the  geodesic  of  the  cone 
(33)  from  P  (  x  ,  y )  =  ( - 1  ,  - 1  )  to  (1,1)  r0  and  0O  in  equation 
(37)  must  be  found.  It  is  argued  that  the  geodesic  from 
P ( x , y )  =  ( -  1 , 0 )  to  (  +  1,0)  has  the  same  length.  For  this  case 
r  =1.6414, 

0  =  11/2. 

The  length  for  this  geodesic  is  s=3.901. 


